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A BST FACT 



Generalized Lancfcester- t ype differential equations are 
used to study combat processes. This system of non-linear 
equations has multiple equilibrium solutions which can be 
determined by a numerical technique called the Continuation 
Method. Useful properties pertaining to neighborhood 
stability are derived by considering the lowest-dimensional 
(1*1) problem. A new set of parameters based on the system 
asymptotes is defined and used to characterize stability. 
System dynamics are investigated using phase trajectories 
which are found to depend on the domains of attraction and 
stabilities of surrounding equilibria. The effect of varying 
initial force levels (X,Y) is studied by calculating an 
objective function which is the difference of the losses at 
the end of a multistage battle simulation. Based on the 
minimax theorem, a set of mixed strategies for (X,Y) can be 
found. For highly unstable warfare with large war resources, 
instability can be used to influence battle outcome. 
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I. IN TRODOCTION 

Since World War II, comtat modeling, simulation and 
analysis have been the subjects of considerable research. 
The objectives of this research are to support defense deci- 
sion making and doctrinal developments during peace and war 
time. During peacetime defense-planners are primarily 
concerned with weapon procurement, development, acquisition, 
organisation and structuring. During war time it is 
believed that a better understanding of the quantitative 
aspects of attrition can help commanders make better command 
and control decisions. 

Ccmbat processes involve complicated interactions 
between opposing forces. These interactions are often influ- 
enced by many external factors such as environment, troop 
quality and tactics. There are different types of ccmbat 
models such as war games, simulations and analytical models. 
A fundamental requirement for a good model is that it must 
be of a fairly high degree of operational realism, since 
otherwise they would not be credible to military planners. 
Cn the other hand, excessively complicated models can make 
the mathematics too difficult tc handle. 

In this thesis, a generalised Lanchester [Ref. 1] model 
which contains area-fire, aimed-fire, self-attrition and 
replenishment coefficients is used. It consists of a system 
of 2N bilinear equations and belongs to the general category 
of analytical models. The model is rich enough to treat 
modern combined-arms operations involving heterogeneous 
forces. It is also possible to extend the model to analyse 
operations cn two or more fronts. 

Among the many important issues that could be analysed 
using this model, the problem of optimum force distribution 



had been studied by fiozencraft and Moose (1983). In their 
paper [Ref. 2], an objective function was chosen as the 
difference of the aggregate attrition rates. It was shewn 
that the optimization problem is mathematically equivalent 
to a matrix game. Hence, the model has a saddle-point solu- 
tion with corresponding optimum force distribution vectors 
x and y for Blue and Crange forces respectively. 

In addition, the neighborhood stability of the model at 
the operating point (x* and y *) was also investigated. By 
defining two parameters, K1 and K2 which are obtained by 
considering small perturbations around the operating points, 
a great deal could he learned about stability. 

Motivated by these results, much of the work done during 
the initial part of this thesis was directed at studying the 
effect of stability cn battle outcome. The ultimate guestion 
is, hew do we exploit the knowledge of stability of an oper- 
ating point to influence battle outcome? Before this gues- 
tion can be answered, it appears that there is a need for a 
better understanding of the equilibrium points. Chapter III 
is devoted to finding and understanding the equilibrium 
solutions and their stability behavior. Like many other 
nonlinear system of equations, the Lanchester's model 
adopted here has multiple equilibria. Stability analysis 
[fief. 3 ] of a non-linear system is usually done by methods 
which do net require prior knowledge of the equilibrium 
solutions. One example of such a method is the Liapunov 
method [Ref. 4], If, by some realizable means, the equilib- 
rium solutions can be found explicitly then there is no need 
to rely on these indirect methods which are often difficult 
to implement. 

One of the reasons for resorting to the Liapunov method 
is the difficulty in obtaining equilibrium solutions of a 
non-linear system. Many numerical methods are unsuitable for 
reasons such as difficulty in obtaining good initial 
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guesses, non-convergence, ill-conditioning and so forth. 
Fortunately, a powerful numerical technique called the 
Continuation Method can be applied for our purpose. This 
method not only finds all the solutions (i.e. it is exhaus- 
tive), it dees not even require initial guesses. 

In order to gain a firm grasp on the dynamics of the 
system surrounding the equilibria, it is helpful to tempo- 
rarily focus attention on the homogeneous (1*1) system. In 
spite of its simplicity, the 1*1 system is not devoid of the 
essential characteristics of the N* N system. In fact, the 
1*1 model is sufficiently sophisticated for certain analyses 
in which the opposing forces can be assumed to be homoge- 
neous. As we proceed through Chapter IV, it will become 
clear that much insight into the stability and system 
dynamics could be gained by merely considering the 1*1 
system. Part of the chapter is devoted to the derivations 
and interpretations of the relations between system asymp- 
totes, locations of equilibrium points and stability. The 
dynamics of the system are studied using the idea of phase 
trajectories. These trajectories represent changes of force 
levels with time and they will be shown to depend not only 
on the stabilities of equilibrium points but also on the 
domains of attraction. 

Chapter V concentrates on battle outcome which is ere of 
the main issues facing a commander. It encompasses many 
issues such as, (1) Who will win and by what margin? (2) 
What is the length of battle? (3) How do initial deployments 
affect battle outcome? (4) Which parameters af feet battle 
outcome most? But we will only address the two following 
subjects : 



(a) The effect 

(b) The effect 
levels . 



of stability on battle outcome; 
of varying X and Y , 



14 



the initial force 



The basic approach is to define a multistage battle with 
a predetermined condition for termination. The resultant 
payoff matrix can then be used to obtain the optimum set of 
mixed strategies. An example, which employs KOREAN WAR data, 
is presented for the purpose of illustrations and 
discussions. 

The essence of tie findings are: 

1. Unstable operating conditions can be exploited to 
influence battle outcome, especially when total war 
resources are large. The effect on battle outcome is 
mere pronounced for highly unstable warfare; 

2. Initial force deployment can be optimized in accor- 
dance with a set of mixed strategies. 

We conclude this introduction by stating two of the 
outstanding issues. The first question is the extent to 
which one can replace the N*N problem by the 1*1 problem. 
The motivation to find an equivalent 1*1 system stems from 
(1) our better understanding of the 1*1 system, (2) ease of 
presenting and visualizing two-dimensional pictures, and (3) 
savings in computational effort. 

The second question concerns replenishment rates. In 
this thesis, the replenishment terms used in the model have 
teen constant. It is therefore reasonable to ask, how to 
modify replenishment terms to reflect a higher degree of 
operational realism? In other words, are there more suit- 

able time-dependent replenishment rates r(t)? 
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II. IA NCH E STEE 1 5 EQUATION 



A. BACKGROUND 

Combat models have been studied as a form of decision 
aid for defense planring. A wide variety of defense plan- 
ning problems, ranging from force structuring and weapon 
selection to rates of deployment in battles have been anal- 
ysed using combat models. There are many different types of 
models. They can be loosely categorized as either war games, 
simulations or analytical models. Discussions on the 
nature, advantages and shortcomings of each can be found in 
[Ref. 5]. 

Cur attention will be focused on a generalized 
lanchester’s [Ref. 5] model, which is an analytical model. 
It consists basically of a system of ordinary differential 
equations describing the mutual interactions -between 
opposing combat forces. Although earlier works in 

lanchester’s model [Ref. 6] employed only a few terms in the 
equations, modern high speed computers enable more general- 
ised, realistic and responsive versions to be used. 

Consider a battlefield with opposing forces, Elue and 
Crange, denoted by { x } and {y.} respectively. The 

subscripts i, j refer to the type of forces such as 

infantry, tanks, artillery, etc. A generalised version of 
lanchester’s model given by 



*i = ‘ x i u i ' x il>ij y j ' 2>ij y j + r i 



v . = - V . v . - V ■ ) X • c • • - } x . d . . + s 

y ] g y ] 11] Zji 1] 



( eqn 2.1) 
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j 



• • 



1.2 



1 2 
x ? - > 



,1 



where 



u i > 


v . 
3 


a i j ’ 


c . . 
13 


b. . , 
i] 


d . . 
13 


r i ’ 


s . 

3 



self -attrition coefficients 
area-fire attrition coefficients 
aimed-fire attrition coefficients 
replenishment coefficients 



is adopted in this thesis. 

Note that in general I ? J, implying that the force 
compositions may be different for the two sides. It is also 
possible to extend the above formulation to a scenerio 
involving more than one battlefield. 

In the next two sections, the highlights of the work 
done by Wozencraft and Moose (1983) are given. The work done 
in this thesis is a continuation and extention of their 
work. The detailed derivations of the results obtained by 
them can be found in [Eef. 2], and hence are not included 



her e . 



B. OPTIMUM FORCE DISTRIBUTION 

The question of optimum force distribution arises in 
combined-arms operations. The problem is fundamentally 
this: Given aggregated forces X, Y, how should one 

distribute them among the different types x i and , i = 

1,2. ..,1, j = 1,2...,J? Since loss rate is one of the 
fundamental concepts in combat modeling, it is reasonable to 
choose this measure as a starting point. The objective 
function was chosen to be 
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( eG n 2.2) 



M 






s . ) 
1 



Per this choice of M, it was shown that there exists 
optimum force distribution (row and column) vectors x* and 
y* such that for any ether vectors x and y 



A ^ ^ ^ A 

xAy < M < x Ay (eqn 2.3) 



where 

k k * k 

M = x Ay 

A 

A = matrix determined by attrition coefficients and 
the aggregate force levels X and Y 

The resemblance of this result to the Mini max theorem 
[Ref. 7] in matrix games is very striking. Indeed, this 
result holds precisely because n can be written in a form 
mathematically eguivalent tea matrix game. Consequently, 
it is net surprising that one can solve for the optimum 
vector x* and v* by means of a Linear Program. An interac- 
tive program to solve a 2*2 program is given in Appendix A. 

C. NEIGHBORHOOD STABILITY 

Equilibrium conditions can be achieved if the replenish- 
ment rates are chosen to make 

0 

1 , 2 . . . ,1 
1,2. . . , J 
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at x = x * and y = y* . Following the usual approach in the 
analysis of nonlinear system stability, equation 2.1 can 
then be transformed into a system of linear equations. 



<5 x 




<5x 




1 

< PQ 

< < 

1 


• 


= -C 




; C = 




6y 




6y 




D C 



C is called the conflict matrix and its elements are deter- 
mined by the attrition coefficients and the optimum vectors 

A A 

x* and y* . For the system of equations 2.1, A and C are 
diagonal matrices. It was shown that two parameters k k 2 
partially characterize the stability of the system. k 1 and k 2 
turn out to be the column sums of the left and right side of 
the matrix 




A A A /\ A 

Eenoting the elements of the submatrices A, B, C, D, by a i -, 
bjj , Cjj and d 1:L respectively, k x and k 2 can be written as 

A \ A A 

k . = - a . . +/ d . . 

1 li / ^ li 



k 2 = -c 



3 3 lj 



independent of the columns i, j. Furthermore, it was shown 
that the following relation holds 



6X - kjSX = Si - k,SY 
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where 



SX » V\5x. ; 6 Y £ V5v. 

l ^ ^ j 

1 j 

It was found that the equilibrium point (x*, y* ) is 
stable if k-^ and k 2 are negative. If kj^ and k 9 are positive, 
then the system is ’unstable*. * Furthermore, values of k 1 
and k 2 and hence the stability of the operating point was 
found to be affected by the aggregate X and Y. 



l Kore generally, it can be_shown that k-, < A <k 9 , where A 
is the maximum eigenvalue of -C. ± u a u 
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III. MU11 IDIM E NSIONAL (N *N) SYSTEM 
A. MATURE CF N*N PROEIEM 

The interesting results highlighted in the last chapter 
provided motivation tc extend the body of knowledge. A study 
of the effect on stability of battle outcome seems to have 
important potentials for applications. Should a commander 
strive to establish a stable operating point, and if so, 
under what conditions? Also, what is the optimum initial 
level of forces he should deploy and how many should he 
maintain in reserve? To answer these questions, more knowl- 
edge about the nature of these equilibria and their 
stability behavior is required. 

The next section outlines the kind of problems we would 
expect to see and their potential complexity. It is 
followed by. a section on finding the equilibrium solutions. 

1 . Exi ste nce of Multi p le Equi lib ria 



An N*N system is in equilibrium if the replenishment 



rates ^ , s. are such that there is no change in the force 
levels {x i - Yj ~ °) • Ihe system of equations becomes 


0 = -x.u. - x . V a . . y . - 

ii i ^ i j 1 j 


Vb. . y . + r . 
^ ij ; j i 


j 


j 

( eqn 3.1) 


0 = - v . y . - y • Vx • c . . - 

r j i ij 


V x . d • . + s . 
^ i ij J 


i 


i 


i > j = 


1 , 2 , . . . , N 
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where, for simplicity, i 
types cf forces. 

A 2N-tuple vector, 
tion 3.1 is an equilibrium 
systems cf equations, equat 
librium point. Geometrica 
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Figure 3.2 Infinite Number of Equilibria. 

in an N*N problem? The answer to this question is not imme- 
diately obvious just by looking at equation 3.1; however, it 
emerges guite naturally when the Continuation Method is 
considered in Section IIIB. It will be seen then that an 
N*N system has, in general, N k equilibrium points where 



N 




Two exceptions, or degenerate cases, have been 
observed, namely; (1) when some or all of the hy per surfaces 
merge there are an infinite number of equilibrium points, 
(see Figure 3.2), (2) when some or all of the hypersurfaces 

intersect in such a manner that repeated equilibria are 
formed, the number of distinct equilibria is less than N k . 
Figure 3.3 illustrates such a degeneracy. 
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Figure 3.3 Repeated Equilibria. 
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Figure 3.4 Domains of Attraction. 

Domains of attraction are separated by boundaries 
which are invariant curves in 1*1 problems and invariant 
hypersurf aces in N * N problems. A boundary surface may be 
considered as an infinite number of invariant curves placed 
side by side. A boundary curve is the locus of points that 
approach an unstable point from both sides. The boundary 
line can be obtained by backward integration (i.e. using 
negative time in equation 2.1) starting just on either side 
of an unstable point. The rationale behind this method is 
that to approach an unstable equilibrium, a point must 
remain exactly on the boundary. If this is not the case, 
then the point will be attracted into the domains and move 
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toward a stable point or infinity. By performing a backward 
integration, we are actually retracing the path taken by a 
point which previously approached the unstable equilibrium 
point. This method requires knowledge of the unstable equi- 
libria, but this is made feasible because the Continuation 
Methods can be used tc find all equilibrium solutions. 

E. IINDING THE EQUILIBRIUM SOLUTIONS 

Tc obtain a set cf equilibrium solutions, one has to 
solve equation 3.1, which can be written using a more 
compact notation as 



F(z) = 0 



( eqn 5.2) 



where 



F(.) represents the right-hand side of 
equation 3.1 



z = (x,y) 

0 = zero vector 



It is well known that numerical techniques for solving 
nonlinear equations are not always successful. Since equa- 
tion 3.2 describes a bilinear system, one should expect to 
face similiar difficulties when attempting to solve it 
numerically. 

Most numerical methods for root finding generally 
require that a fairly good initial guess (z Q ) be known so 
that some convergent iteration process 

*n+l ' 8( V 
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brings the approximated root closer and closer to z the 
desired equilibrium solution or root. In practice, the 
following difficulties are often encountered : 

(1) The convergence condition of the algorithm 
must be ensured ; 

(2) Finding an initial guess that is sufficiently 
close to the correct solution is difficult, 
especially for higher dimensions; 

(3) Even if a good initial guess has been obtained, 
the numerical process may still be plagued by 
ill-conditioning, saddle points, etc.; 

(4) Not all the solutions are guaranteed to be 
foun d. 

^ • Con tin u a ti on M etho d 

Fortunately, the above problems are avoided if a 
numerical method called the Continuation Method [Ref. 8] is 
used. This technique, which is sometimes called The 
Imbedding Method, has been successfully applied in many 
fields. It introduces an artifical guide which will channel 
the iterates toward a specific solution. Such a guiding 
principle is actually a knowledge of the existence of a 
suitable curve connecting an initial point with the desired 
solution. 

Continuation Method has significant advantages over 
other numerical techniques. Most importantly, a good initial 
guess is not necessary and all the solutions can be 
obtained. 

a. Basic Theory 

Given the problem F(z) = 0 to solve, the first 
step is to embed it into a homotopy or a parameterized set 
of problems, H (z,t) . The requirements on H(z,t) are : 
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(1) 9(2,1) = F-l (z) =0 is the original problem 

(2) H (z ,0) = F q (z) = 0 has a trivial or easily 

computed solution 

For example, a homotopy could be : 

H ( z , t ) - tF ( z) + (l-t)F 0 (z) , t€[0,l] (eqn. 3. 3) 

Using the above parameterization, the simple 
problem of F 0 (z) = 0 is deformed into the desired one, Fj_ (z) 

= 0 . This is done by calculating the solution tc the 

deformed problem at each stage of the deformation. The exis- 
tence of a continuous curve such that H(z(t), t) is a solu- 
tion to H(.,.) = 0 for all t£[0, 1] is assumed. 

b. Implementation 

To actually carry out the above continuation 
process one usually differentiates H(.,.) to form 

H(z(t) , t) = 0 (eqn 3.4) 



Using eguation 3.4, z can be written as a function of z and 
t as given in equation 3.5. The function, h(.,.) is prefer- 
ably a linear function that can be integrated numerically. 

i = h (z , t) (eqn 3 . 5) 

Together with the initial condition z (0) = z Q , 

equation 3.5 is actually an initial value problem which can 
be integrated numerically. The solution at t=1 is then the 
solution to the original problem F(z) = 0. 
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2. Algorithm to Obtai n 2* 2 Equilibrium P r ob i e m 



A 2*2 Lanchester problem is first formulated into a 
Continuation process. It is followed by a discussion on how 
the accuracy of the method can be improved. The last part of 
this subsection includes a note on the number of equilibrium 
points in an N*N problem. 

a. Formulation 

For the 2*2 problem, F (z) =0 is explictly 



- z l (u l +a ll z 3 +a 12 z 4> 


+ 


r i 


" b l 1 Z 3 


b 1 2 Z 4 


" Z 2 ( u 2 + a 21 Z 3 + a 22 z 4' ) 


+ 


r 2 


" b 2 1 z 3 


b 2 2 Z 4 


' Z 3^ U 3 +C 11 Z 1 +C 21 Z 2^ 


+ 


r 3 


- r\ *7 

a ll“l 


d 2 1 z 2 


' z 4 (u 4 + C 12 Z l " C 22 Z 2 ) 


+ 


r 4 


“ d. 7 
12**! 


d 2 2 Z 2 



The homotopy is formed by writing 

H 1 (z) = H 1 (z,0) + r ^ r i' b n z 3" b i2 z 4^ = 0 

H ? (z) = ( z , 0) + t(r 0 -b o .z -b ?9 z.) = 0 

2 " Ll 0 2i * 4 (eqn 3.6) 

H 3 (z) = H 3 (z,0) + t (r 3 -d ll z l' d 21 Z 2' ) = 0 

H 4 (z) = H 4 (z, 0) + t (r 4' d 12 z i' d 22 Z 2^ = 0 

where 



= 


' Z 1 ^ u l +a ll z 3 +a 12 Z 4^ 


= 0 


H 2 (z,0) = 


' z 2 (u 2 +a 21 Z 3 +a 22 Z 4^ 


= 0 


HjCz.O) = 


' Z 3 (U 3 +C 11 Z 1 +C 21 Z 2^ 


= 0 


H 4 (z,0) - 


" Z 4 ^ U 4 +C 12 Z 1 +C 22 Z 2^ 


= 0 
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Next, we differentiate equation 3.6 with respect to t and 
put it in a matrix form 



Az = B (eqn 3.7) 

where 



r- 








V a llV a 12 z 4 


0 


a ll Z l +b ll t 


a 12 Z l +b 12 t 


0 


u 2 +a 21 Z 3 +a 22 Z 4 


a 21 Z 2 +b 21 t 


a 22 Z 2 +b 22 t 


c ir3 +d n t 


C 21 Z 3 +d 21 t 


U 3 +C 11 Z 1 +C 21 Z 


2 0 


C 12 : 4 +d 12 t 


C 22 Z 4 +d 22 t 


0 


u.+c, -,z. +c 70 z 9 
4 12 1 22 2j 



z [ z i» z 2» z 3» z 4^ 



B = 



r l‘ b ll Z 5' b 12 Z 4 

r 2' b 21 Z 3" b 22 Z 4 

r 3' d ll Z l' d 21 Z 2 

r 4' d 12 Z l' d 22 Z 2 



Equation 3.7 can now he integrated numerically using one of 
the readily available integration routines. 

tfe have assumed that the trivial solution to 
H(z, 0) s 0 has been previously found. 

b. Improving Accuracy 

Numerical integration of equation 3.7 inevitatly 
produces seme errors at each iteration. Since the 
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Continuation method relies on following curves to arrive at 
the desired solution, it is esssential that each iterate 
remains close to the actual curve. It is necessary to 
include a way to correct the approximated position fcy means 
cf a corrector step. The combination of integration and 
correction is often called a "predictor-corrector step" 

This process of prediction-correction is shewn 
in Figure 3.5 where each integration error has been exagger- 
ated for illustrative purposes. The algorithm tc be 
presented later employs an IMSL routine called ZSCNT for the 
predictor step. Other forms of curve following routine can 
also be found in the literature, and are briefly mentioned 
in [Eef. 8]. 

c. Trivial Solution 



The trivial system H (z, 0) = 0 was chosen to be 



H 1 


(2,0) 


- Z 1 ( u i + a n z 3 + a i2 z 4^ 


= 0 




H 2 


(z , 0) 


' z 2 (u 2 +a 21 Z 3 +a 22 Z 4^ 


= 0 


(eqn 3.8) 


H 3 


(2,0) 


‘ Z 3 ( - U 3 + C 11 Z 1 + C 21 Z 2^ 


= 0 




H 4 


(2,0) 


~ z 4 ( U 4 + C 12 Z 1 + C 22 Z 2' ) 


= 0 






In 


non- d egenera te cases. 


there 


are six scluticns 



corresponding to equation 3.8. The result is derived in 
Appendix B which also deduces the number of trivial solu- 
tions for an N*N problem to be 




(eqn 3.9) 



of obtaining the trivial solutions is given in 
Using a combinatorial identity, N, can be 

K 



3 1 



The method 
Appendix 3. 
written as 




Figure 3.5 Integration Path using Predictor-Corrector. 

/2N 

N * = (n 

Each continuation process starts from a trivial 
solution 2 Q and follows a specific curve until it reaches 
the equilibrium point. Consequently, the number of equilib- 
rium points will also be N k , As mentioned in section IIIA, 
the two exceptions are situations involving inf initely-many 
and repeated equilibria. Situations involving degeneracy 
are discussed in Appendix B. 

d. Algorithm 

(1) Singularity Treatment . In Continuation 
Method algorithms [Ref. 8], it is sometimes necessary to 
give special treatment to cases in which the curves being 
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followed by the integration routine pass through a singu- 
larity. Experimentally, it had been observed that in our 
problem, the singularity took on the form shown in Figure 
3.6. Corrective measures were necessary to ensure that upon 
crossing the singularity, the large magnitude was preserved 
but the sign was changed; otherwise the curve might termi- 
nate at an equilibrium point which was not the intended one. 

In the algorithm, the presence of the 
singularity is detected by monitoring the rate of change of 
the individual component z ± . Once identified, this fast- 
changing and large-magnitude component (z p ) is monitored at 
each step t where 0 = t Q <-..<... t k < t k <...< t 0nd = 1- 
When z F is found not to cross the singularity and end up at 
approximately -z p , the algorithm attempts to correct this 
irregularity by artificially making z p = -z before the next 
predictor step commences. 

(2) Flowchart. The flowchart for the algo- 
rithm is given in figure 3.7. Only the major steps have been 
shown. The program listing is given in Appendix C. 

3 . Example and Ee sults 



Consider as an example a 2*2 problem with the 
following attrition coefficients 
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Figure 3.6 Curve Passing through Singularity. 

The trivial solutions are first computed and serve 
as one of the inputs to the program. The program obtains the 
values of z (t) and plots each component ( (t) ) versus t. 

In Figure 3.8, the plots for t close to zero show one set 
of curves for z i (t) starting from their respective trivial 
solutions. The curves of zp(t) versus t for all the six 
sets cf equilibrium solutions are shown in Figure 3.9. 

A few interesting features of the continuation 
process are worth noting. For example 

• Each trivial solution leads to different equilib- 
riuii solution and the integration path is different 
for each component. 

• All the curves are smooth ; one of the four curves 

may pass through a singularity. ( see Figure 3.9 (c) 

and (d) ) . 

Table I summarizes the computed equilibrium solu- 
tions. They are tabulated in the same order as the plots in 
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(Start) 




Legend 



TS = partition 0<t<l ; A = integration step; 

_N = number of equations ; omax = a constant; 

z - trivial solution, 
o 



F igure 3 . 7 



Flowchart for Algorithm to find 2 : ' ; 2 Equilibria 
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o.l 



Z 3 (f) 




> f 



Figure 3-8 z ^ (t) Versus t for Values of t Close to Zero. 

Figure 3.9. To estimate the accuracies of the results, we 
defined error as 

V ~ 2 “* 2 -2 ~ 2 

F-, + F ? + F, + F, 

1 L o 4 



where 

“ F^(z) » z is the comput ed equilibrium 
solution to F(z) = 0 

Decreasing the integration step size in the 
predictor and corrector routines may reduce the errors by a 
small amount; but the increase in computational effort may 
not be justifiable. Conversely, it may be desirable to cut 
down computing time. Currently, the algorithm performs one 
corrector step for each predictor step. If two or more 
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Figure 3 






Plots of z i (t) Versus t during Continuation Process 
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Figure 3.9 (contd.) 
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TABLE I 



Computed Equilibria for X. = 1.0, Y = 4 . 0 



Trivial Solution 


Computed EauilLtrium Solution, z 


Error 


(0, 0, 
0, 0) 


(1.7755, 1.2906, 
1.0533, 0.3512 


0 . 2*10" 4 


(0 , -0.33, 
0 , -0.33) 


(34.8989, -64.2951, 
-0.1350 , -0 . 2776 ) 


0 .6*10"° 


(0, -0.091, 
-0.5, 0 ) 


( 0.615 4 , 0 . 394 5 , 
3 .0769 , 0 . 9 2 31 ) 


0 . 1*10”° 


(-0.2, 0, 
0, -1.0) 


(0.46- : 10" 4 , -0.39 30 , 
-16.6400, 52.9920 ) 


0 . 22* 10" 5 


(-0.083, 0, 
-0.3, 0) 


(-0.135 7 , 0 . 0 32 5 , 
-44.5297, 28.6147) 


0 . 56- ; 10 _6 


(-0.42, 0.37, 
0.25, -0.17) 


(-45.2284, 39.1604, 

-0.3497, -0.04485) 


0 . 2 3-10 -6 



predictor steps are done for each corrector step, seme 
computational effort 2 can be saved. 



2 Saving in computational effort will be more significant 
when solving higher dimensional systems. 
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IV. PBCEE RTIE S OF THE JjM SYSTEM 



The 1*1 problem is the simplest case in our model. It is 
nevertheless important for us to investigate and understand 
its properties. Despite its relative simplicity, it is by no 
means uninteresting. There exists many situations which can 
he realistically and easily modeled by the 1*1 system. For 
example, when the opposing forces can be considered as homo- 
geneous, it is convenient to use the 1*1 model for analysis. 
It is also useful for the analyses at the strategic level 
when the forces and parameters can be aggregated. In many 
instances, it seems to provide insight on how to approach 
the N*N problem, which is much more difficult to visualize. 
In fact, as the understanding of the 1*1 system increases, 
there is a strong urge to try to represent the N*N problem 
by an equivalent 1*1 problem. The equivalent representation 
is not only attractive in terms of its simplicity but also 
its economy in computational efforts. 

The next section will focus on the relation between 
system asymptotes and stability of the equilibria. By formu- 
lating the problem quantitatively, we are able to arrive at 
some useful properties. In Section IVB, the system dynamics 
i.e. the changes in the force levels are analysed by consid- 
ering the phase trajectories. 



A. SYSTEM ASYMPTOTES AND EQUIIIBRIUM POINTS 
For the 1*1 problem, the system reduces to 



x = -x(u + ay) + r - by 
y = -y(v + cx) + s - dx 



(eqn 4.1) 
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An equilibrium condition exists if r and s are chosen such 
that x and y are both zero. In general, there will be two 
equilibrium points corresponding to two locations where the 
two hyperbolas intersect. The hyperbolas are described by 



x 



r - by 
u + ay 



s - dx 
v + cx 



(eqn 4.2) 



From equation 4. 2, one can easily deduce the fouc asymptotes 
(two vertical and two horizontal) associated with the hyper- 
bolas. Figure 4.1 shows a typical set of four asymptotes. 
They always cross in the third guadrant of the x- y plane and 
do not depend on the replenishment coefficients. The rela- 
tive displacement s between the two horizontal (and also 
vertical) asymptotes depend only on the ratios of attrition 
coefficients and not on the coefficients themselves. It 
turns out that these properties of the system asymptotes 
help to simplify the analysis considerably . 



1 . Sta bil ity Criteria 

Considering small perturbations about an equilibrium 
(x e ,y e ) and linearizing the equations, we have 




(u + ay e ) (b + ax e ) 




6x 


(d + cy e ) (v + cx e )_ 




6y 



The characteristic polynomial is simply 




Figure 4. 1 System Asymptotes. 



D($) = Det [si - C] (eqn 4 - 3 ) 

where 

I = identity matrix 
C = the 2*2 matrix in equation 4.3 

Hence , 

D (s) = [s + (u + ay e )J^s + O+cx^ - [(b + ax e ) (d + cy e )] 

= s 2 + £(u+ay e ) + (v + cx e ) s + j\u + ay e ) (v + cx e ) 

- (b + ax e ) (d + cy e )J 

Ihe conditions for (x e , y e ) to be a stable equilibrium, 
i.e. for the roots of D(s) to be in the Left Half Plane 
(LHP) are given by 



4 2 



(eqn 4.4) 



(u + ay e ) + (v + cx g ) > 0 

(u + ay e ) (v + c x e ) - (b + ax g ) (d + cy Q ) > 0 

2 • Stability and A sym pt otes 

Five different ways in which the hyparbolas can 
intercept have been identified and their stabilities 
accounted for. These five cases are shown in Figure 4.2 and 
each case will be elaborated upon subse quen tly. 

a. Definitions and Formulation 

One of the most intriguing facets of the 1*1 
problem is the connection between the asymptotes and the 
stability cf the resulting equilibria. We begin the quanti- 
tative treatment by first defining the following ratios: 



n = H r) 4 d 

la* 2 c 




The four asymptotes are x = -p 1 , x = -u 2 , y =-r i 1 and y = 

-n 2 • If we let the first equilibrium point be (x el , y el ) 
and substitute the corresponding r and s into the equation 
4.2, we have 

V x ' x el ) + (xy ' X el y el ) + “l^el 5 = 0 

(eqn 4.5) 

n 2 (x-x el ) + (xy-x el y el ) * V^el 5 = 0 
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case (e) 



unstable 

stable 

neutrally stable 



Figure 4.2 Types of Equilibria 
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Next, the distances between the asymptotes are defined as 



£ 

x 



A 



e 

y 



A 




P 



1 




n 



1 



It is not difficult to see that e x and e y will 
decide where the hyperbolas intersect. For instance, when. e x 
> 0 and £ y > 0, there may be two equilibria in the first 
quadrant 3 (See case (b) of Figure 4.2). In general, the 
second equilibrium point ( x e2 / Y e2 ) can be b° un d by elimi- 
nating y or x from equation 4.5 and comparing coefficients 
with (y - y el ) (y - y e2 ) and (x - x el ) (x - x &2 ) . The final 
expressions are 



x , = — (y . + n_) - u 
e2 e w el 1 

e 

y - = — (x + P . ) - n 
7 e 2 e v el 1 



(eqn 4.6) 



For constant x el , x e2 / y el # an ^ Y e 2 ' equation 
4.6 can be written to represent two straight lines in £ x , 
e y The equations of these two lines are 



e 

y 



e 



y 



£ 

X 



£ 

x 



cy e z + V 

Uei ♦ V 

!VV 

(x e2 + 



(eqn 4.7) 



3 In our context, the quadrants are defined by the asymp- 
totes and net by the x, y axes. 



b. Types of Equilibria and their Stability 

To derive the different types of equilibria and 
their associated stabilities, we make a transition from the 
x, y plane into the e x , e y plane. Briefly, the basic 
approach is to fix ore equilibrium point (x , , y ) on the 
first guadrant hyperbolas and consider the regions in the 
£ x # Cy Pi ane when we have the other point in various places 
of the x, y plane. The other essential step is to express 
the stability criteria (equation 4.4) in terms of e x , e , 
P 1# n 1 , x el ' y e i • A summary of the results which are derived 
in Appendix D is given below ; 

(1) When both equilibrium points are on the first 
guadrant hyperbolas ( case (b) in Figure 4.2 ), one 
will be stable and the other unstable; 

(2) When one equilibrium point is on the first 
guadrant and the other on the third, both can be 
unstable or one will be stable and the other 
unstable ( case (a) in Figure 4.2 ) ; 

(3) When both equilibrium points are on the third 
guadrant hyperbolas, both are unstable ( case (c) in 
Figure 4.2 ) ; 

(4) When there are infinite number of equilibria as 
in case (d) in Figure 4.2, £ x = e y = 0 and the two 
sets of hyperbolas merge. Equilibria lying on the 
first quadrant hyperbola are neutrally stable (one 
eigenvalue equals zero) and those on the ether 
hyperbola are unstable; 

(5) When there are repeated equilibria as in case 
(e) in Figure 4.2, they are neutrally stable if the 
hyperbolas touch in the first guadrant ; otherwise 
they are unstable. 
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Host of the above results are embedded within 
Figure 4.3 which is reproduced from Appendix D for conven- 
ience. Evidently/ both the coordinates of the eguilibrium 
points (x e , y e ) and the location in the e x , e y plane deter- 
mine the stabilities. The e x , e plane has been subdivided 
into a few regions each with distinct stability 
char act er is tics. 





^ x el ’ y eP 
unstable; 
(x e2 , y e2 ) 

stable 




as above 
both 

unstable 




(x el’ y el 
stable; 

(x e2> >eZ 
unstable 



) 

) 




as above 



c y e x' 







) 

) 



Figure 4.3 The e , e Plane. 

x y 

The case of infinitely many eguilibria corre- 
sponds to the origin cf e x , £y plane ( V 1 = b 2 0 = n 2 ) . The 

only way for two sets of hyperbolas to merge is for their 
respective asymptotes to merge. This case is a degenerate 
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instance of repeated equilibria ( case (e) in Figure 4.2 ), 
which is shown in Appendix D to correspond to operating 
points on the line e ^ = e ^ (Y + n J / ( X + p 1 ) as illustrated 
in Figure 4.3. 

As a corollary, we note that there cannot be two 
stable equilibrium pcints in the 1*1 problem. This deduction 
can be made by referring to Figure 4.3. There is no region 
in the £ x , £ y plane which allows for this case. At most, 
there can be two neutral ly stable equilibria which are 
repeated. Numerous attempts have been made to obtain two 
stable equilibria in the 2*2 problem, but in vain. Whether 
it is also true for 2*2 or higher dimensional problems that 
only cne equilibrium may be stable is still a matter of 
conjecture. 

In Appendix E, the relations between the regions 
on the e x , e y plane and their associated stabilities are 
verified. Some representative points on the e x , e y plane 
are chosen and their stabilities checked. 

E. SYSTEM DYNAMICS 

The dynamics of a 1*1 system are characterised by its 
phase trajectories, which are curves on the x-y plane 
describing the history of the system as the time, t, 
changes. These trajectories can be conveniently obtained by 
integrating equation 4.1 numerically. 

Needless to say, being able to predict the trajectories 
is important, for it means that we know how our model of a 
tattle progresses. Cnee the factors influencing the course 
of a battle are known, appropriate command decisions can be 
introduced to ensure favorable battle outcome. In Chapter V, 
we will see how many cf the results obtained in this section 
can be used to rationalize and predict battle outcome. 
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Some typical trajectories corr espondin 9 to the different 
types of equilibria are described in the next subsection. 
Besides the stability which influences trajectories, it was 
briefly mentioned in Chapter III that domains of attraction 
also affect the trajectories. In the subsection that 
follows, we will show specific examples of the way to deter- 
mine the domains by finding their exact boundaries. 

1. Ira iect ori es 

Two methods of establishing the trajectories from a 
given initial condition will be described. The brute-force 
method which has been mentioned uses numerical integration. 
The ether method which often provides better insight, is 
more graphical. The graphical method is based on a few very 
simple rules to predict the gross behavior of a trajectory. 
Some of these rules are listed below : 

(1) A stable point "attracts"; unstable point 

"repels" ; 

(2) Points on either side of a boundary move into 

their respective domains; 

(3) For large (x, y) , trajectories are governed by 

the Lanchester "linear law"; 

(4) Points near the hyperbolas can be easily 

• • 

analyzed by noting the signs of x and y. 

As an example of using the graphical method to 
determine trajectories, consider a region around an unstable 
equilibrium point on the first quadrant hyperbola. The whole 
picture of the phase trajectories (sometimes called phase- 
plane portrait [Ref. 9] ) can be put together in a logical 
fashion by using those simple rules. Since this equilibrium 
point is unstable, trajectories will be expected to diverge 
from it. As an unstable equilibrium point, it will have a 
boundary line passing through it. Initial conditions start 



49 



from each side give rise to different trajectories. Next, we 
determine the signs cf x , y on both sides of each hyperbola 
as indicated in Figure 4.4 where only one intersection is 
shown . 



y 




Figure 4,4 Analytical method of predicting trajectories. 



Note how predictable these trajectories are. If, 
for some reasons, the exact trajectories are required, we 
can resort to the brute-force method. The methods are obvi- 
ously complementary in nature. The advantages of the brute- 
force method are accuracy and simplicity. In Figure 4.5, a 
typical computer plot consisting of ten trajectories is 
shown. The program which produces the plot is included in 
Appendix F. 

Referring to Figure 4.5, the trajectories cross the 
hyperbolas and move asymptotically along a common curve 
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y 




Figure 4.5 Computer Plot of Trajectories. 

lying between the hyperbolas. This same property is exhib- 
ited by other cases. Even the special case with no hyper- 
bolic intersection has been found to behave similarly as can 
seen in Figure 4.6. 

Our ability to determine the trajectories and 
present them vividly is partly due to fact that two- 
dimensional pictures can be easily drawn and visualized. For 
dimensions higher than the third, it is impossible to visu- 
alize trajectories; however, the notion of trajectories can 
be conceptually extended to n-dimensional space. Thus, it 
seems likely that in the higher dimensional systems, trajec- 
tories cross hypersurfaces and move along a common asymp- 
totic curve analogous to that in the 1*1 system. Further 
studies are required before this behavior can be confirmed. 
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Figure 4.6 Trajectories when Hyperbolas do not Intersect. 

2 • Bo undaries of Doma i ns of Attraction 

In Chapter III, the idea of the domains of attrac- 
tion was briefly discussed. In an n-dimen sional space, such 
a domain is a region or volume in which all initial points 
come under similiar influence. When domains exist, there 
will be boundary surfaces which can be thought cf as 
collections of invariant curves passing through unstable 
equilibria. 

For a 1*1 problem, domains and boundaries are net at 
all abstract. In the last subsection, they have been shown 
to affect trajectories. Recall that in Chapter III, we 
mentioned a simple and yet effective way of finding the 
boundary curves and establishing the domains in the x-y 
plane. Examples on the use of backward integration to obtain 
boundary curves are now presented. 
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rting from an unstable point, we apply small 
n both directions perpendicular to an eigen- 
ted with the most positive eigenvalue and 
ard in time (in the computer program, this is 
employing negative time steps for integra- 
sult is a smooth, invariant curve which is 
undary or the so-called separatrix like the 
gure 4.7. 




Figure 4.7 Boundary Curve through an Unstable Point. 

To verify that the curve is indeed the boundary, 
two initial points are chosen just off the curve (e.g A, B 
in Figure 4.7) . if we forward integrate from these two 
points, they move into different domains as indicated in the 
same diagram. Appendix G contains a Fortran program that 
does the backward integration and plots the boundary curve. 
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3. Bou nda ry Curve between Two Hyper bo las 

Boundary curves do not necessarily pass through 
unstable points. Backward integration methods can also be 
used if a boundary exists but there is no unstable equilib- 
rium point to serve as the starting point of integration. 
This is best illustrated by considering the case of both 
equilibria cn the first quadrant hyperbolas. In this case, 
there is no equilibrium point in the third quadrant; never- 
theless a boundary does exist between the third-quadrant 
hyperbolas. The existence of the boundary is visible by 
simply considering the signs of x and y on both sides of the 
hyperbolas. In figure 4.8, the signs of x and y and also the 
directions of some typical trajectories are depicted. 




Figure 4.8 Existence of Boundary Between Two Hyperbolas. 

To obtain the exact boundary, choose a point close 
to a hyperbola and cn lower part of the hyperbolas (e.g. 
point P or Q in Figure 4. 8) and integrate backward. The 
result is a boundary curve as shown in Figure 4.9. 
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Figure 4.9 Exact Boundary Curve Between Two Hyperbolas. 

4 . Sum mar y of t he 1*1 P robl em 

We have seen the close relation between system 
asymptotes and stabilities. Through the use of newly defined 
variables e x and e , the stability of different types of 
equilibria has been derived. Five cases have been identi- 
fied, and they correspond to the types of inter sections on 
the x-y plane. For example, if both the equilibria are found 
on the third guadrant hyperbolas, then we know that they 
will be unstable. 

Two methods of establishing the trajectories have 
been described in this chapter. These two methods complement 
each other and the choice depends on our requirements. Ihe 
dynamics of the system are characterized by the trajecto- 
ries, which as we have seen are very predictable. These 
trajectories are influenced by the stabilities of equilibria 
and domains of attraction which are separated by boundary 
curves. A simple way of plotting the boundary curves has 
also been presented along with specific examples. 
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The results derived in this chapter will be applied 
in the next chapter. The knowledge of the system dynamics 
and how they are affected by stability and other parameters 
will enable us to analyze changes in force levels as the 
battle progresses. 



\ 
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V. STRATEGY FOE INITIAL FORCE COMMITMENT 

Id the last two chapters, emphasis has been placed on 
establishing the mathematical framework of the system 
dynamics and stability. In this chapter, we examine some 
model operational problems that are related to stability and 
dynamic considerations. 

One of the major command decisions that has to be made 
during a build-up period of a war pertains to initial force 
commitment. A good strategy calls for a balance between 
initial deployment and reserves. In practice, a multitude of 
factors have to be considered before deciding on a partic- 
ular commitment. The approach in this chapter provides us 
with a set of mixed strategies but does not consider intan- 
giable factors like world politics, national economy, 
survival factor and so on. 

Stability has been shown to effect trajectories which in 
turn effect battle outcome. Recall from Chapter IV that 
there are some trajectories which represent speedy and 
complete annihilation of one force; hence it seems reason- 
able that the side that is tipped to win the battle will 
want to operate on an unstable trajectory. But to what 
extent can one exploit the stability behavior of the system 
to influence battle outcome? Obviously there will be prac- 
tical limitations; an important one of these is total avail- 
able resources. 

A. PBOBLEH STATEMENT AND APPROACH 

The problem statement is as follows : 

Given total defense resources Q x , Q y for x and y 
respec tivel y, what is the optimum set of strategies for 
initial force commitment, X and Y? 
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W € begin by treating this as a 1 *1 problem at the stra- 
tegic level. The dynamics of the problem are thus governed 
by eguation 4.1. Both sides are assumed to operate initially 
at equilibrium with constant replenishment rates given by 



r = X(u + aY) + bY 
s = y (v + cX) + dX 



(eqn 5.1) 



Since both sides have limited 
the replenishment rates versus 
Figure 5.1. where Q = rT and Q 

3 x x y 



defense resources Q , Q.., 

x y 

time may be as shown in 

= sT . 

y 



replenishment 




Figure 5. 1 Replenishment Versus Time. 

The next step is to select some suitable form of payoff 
function which is to be optimized for a certain choice of X 
and Y. The payoff function (from X to Y) has been chosen to 
be 

A (X , Y) = L y - L x 

where L x , I = Total losses for x, y at battle termination 
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During this stage, the dynamics of the system is 
given by the familiar 1*1 system 



x = -x(u + ay) - by + r 
y = -y(v + cx) - dx + s 



(eqn 5.2) 



Khen this 1*| system is integrated, just as in 
Chapter IV, the resulting trajectories behave similiarly. 
However, there is a major difference. Now, we no longer 
have unlimited defense resources, and this stage will not 
last forever. It implies that, unless Q x or Q y is extremely 
large 4 , trajectories which reflect quick annihilation of 
one of the forces are rare. In general, T x and T y are given 
by 



T 



x 



T 



y 




If one of the force levels drops to less than ten 
percent of Q x or Q y , the battle is arbitrarily considered 
over and the losses are calculated as in Figure 5.2. Ihe 
finish time (FINTIM) is simply t, the time when x < 0. 1Q x or 
Y < 0.lQ y . 

2. Stage 2 

Since either x or y can run out of reserves first, 
the dynamics of stage 2 are governed by either equation 5.3 
cr 5.4 respectively. 



4 Q or Q may be very large if x or y is backed by a 
superpower who is fully committed to provide military aid. 
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Figure 5.2 Losses at Stage 1. 



x = -x(u + ay) - by 
y = -y(v + cx) dx + s 



(eqn 5.3) 



x = -x(u + ay) - by + r 
y = -y (v + cx) - dx 



(eqn 5.4) 



Unless the battle ends earlier, this period will 
last for T 2 which is given by 



T 0 = Max IT , T | - T, 

2 * x y » 1 



During this period, the trajectory will be different from 
that in stage 1. This is because when r = 0 or s = 0, cne of 
the hyperbolas is shifted so as to cross the origin and we 
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have different equilibrium points. The trajectory will now 
be influenced by the new equilibrium point. This is illus- 
trated in Figure 5.3. 

Calculations of the losses are more complex than in 
stage 1 since there are now two cases to deal with i.e. r = 
0 or s = 0. The procedure is shown in Figure 5.4. 

3 . Stage 3 

If the battle enters stage 3 without either x < 
0. 1 Q x or y < 0. IQ then the dynamics will be dominated by 
attritions since r = s = 0. Equation 5.5 is now used for 
integration . 



x = -x(u + ay) - by 
y = -y(v + cx) - dx 



(eqn 5.5) 



Again, the trajectory will have to change because 
now both hyperbolas pass through the origin. This is illus- 
trated in Figure 5.5 where we show how the intersection at 
stage 2 has changed. Losses and FINTIM are calculated in 
accordance with the procedure in Figure 5.6. 



C. MIXED STRATEGIES 

The range 0 to Q 5 for both X and Y can be subdivided 
into id force levels. There are m*m pairs of X and Y and 

corresponding number of payoffs, A(X,Y). We thus have an m*m 
payoff matrix having elements A(X,Y) . Figure 5.7 gives a 
pictorial representation of this two-person game. 



5 In the actual program. one may wish to restrict the 
range of X and Y to interval (0.2Q - 0.75Q) to reflect prac- 
tical limitations in initial force deployment. 
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Figure 5.3 Trajectories and Hyperbolic Intersection During Stage 




Figure 5.4 losses at Stage 2. 

In the last section, the procedure for computing A (X,Y) 
has been described. A simple program can be written to 
compute each element of the payoff matrix. One such program 
is given in Appendix H. 
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Figure 5.5 Trajectories and Hyperbolic Intersection During Stage 



integrate eqn 5 . 5 




Figure 5.6 losses at Stage 3. 

There are a few ways of presenting and interpreting the 
payoff matrix. A normal practice is to present it in tabular 
form and consider only pure strategy. Alternatively, a plot 
of A (X, Y) as a function of X and Y could be obtained. When 
using pure strategies, it has been observed that the game 
does not always have a saddle point [Ref. 10] and it would 
be better to use mixed strategies. 

In mixed strategies, x and y may play all their strat- 
egies in accordance with a certain set of probabilities. 
Although in our situation, x and y can only play once, the 
same concept of mixed strategies is still useful. If we let 

p and g be the probabilities by which x and y select their 

i j 

ith and jth pure strategies respectively, then 
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Figure 5.7 Payoff Matrix 
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In addition the (i,j)th entry of the payoff matrix be 
denoted by aij, the probabilities can be represented by the 
matrix below 
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Ihe optimal mixed strategy is based on the mir.imax 
criterion- Mathematically, x and y select p i and g which 
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will yield U and V as given equation 5.6 and 
respectively. 



max 

p i 
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m 
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£ a ij q j 
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j-i 
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II 

h-* 


j=l J 



equation 5.7 



(eqn 5.6) 



(eqn 5.7) 



Appendix I describes how the problem of solving for the 
optimal values of p ^ and q can be put into linear program- 
ming form. The program given in Appendix H also computes 
this optimum set of solution in addition to obtaining the 
payoff matrix. 

The concept of mixed strategies is quite intuitive if a 
game is to be played repeatedly. But since we are using it 
to provide us with an optimum set of probabilities of 
selecting the pure strategies, some interpretation is 
required. Although the optimum mixed strategies have been 
obtained, a pure strategy still has to be selected and used. 
However, it is important that the selection process should 
be random 6 according to the optimized probabilities 
obtained. 

One simple but valid statistical procedure [Ref. 11] to 
select a pure strategy from a set of mixed strategies is to 
first plot the probability distribution function. A random 
number generator is then used to generate a number between 
zero and one. The corresponding value of the strategy could 
then be selected. This procedure is shown in Figure 5.8. 



6 The selection 
opponent can select 



process must be random otherwise 
a strategy to improve his outcome. 



the 
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Figure 5-8 Obtaining Pure Strategy from Mixed Strategies. 

D- EXAMPLE USING KOREAN SAB DATA 

One cf the main objectives of using actual historical 
data in a model is for validation. It is important that the 
results obtained using the model should at least be consis- 
tent with actual events. The Korean War has been chosen 
because there was a clear-cut victor during the initial 
phase of the war. We consider the period when only North 
Korea and Republic of Korea (South Korea) were involved. 

Before the entire simulation can be carried out, the 
actual force strategies, fighting ability, weapon state, 
etc, have to be transformed into familiar quantities and 
parameters such as Q x , Q y , X, Y, a, b, c, d,..., and so on. 
This transformation, together with some background data on 
the Korean War are given in Appendix J. 

1 • Res ults and Dis cuss ions 

First we examine the resultant trajectories during 
the three stages of the battle which are shown in Figure 
5-9. The simulation uses the X and Y which correspond to 
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the actual initial deployment by both North and South Korea 
respectively. Clearly we see that the victor is x, as it 
was in history. The result of the simulation also shows the 
three stages explained in the last section. Note that the 
trajectories for the first and second stages are curtailed 
because both sides run out of war reserves. The implication 
is that in practice, the kind of trajectories leading to 
large and rapid changes in force levels are rather rare. 

However, the effect of instability on battle outcome 
is borne out by experimenting with the directions of pertur- 
bations. Consider the case in which x (North Korea) fixes 
the initial force and y (South Korea) varying the initial 
force levels around the equilibrium point. In Figure 5.9, 
these perturbed points are denoted by points A to D spanning 
across the boundary separating the domains of attraction. 
From our understanding of the stability and system dynamics 
each perturbation will give rise to different trajectory and 
payoff at the end of the simulation. Clearly, y will want to 
operate at the perturbed points A or B rather C or D since 
the former will result in the trajectory for stage one to be 
in a decreasing x direction. Table II shows the variation in 
the payoff as the perturbation point changes. When the 
perturbed points are at A or E, the payoffs to x are less 
then those for points C or D. Thus we have seen how an 
unstable system can be used to inflict heavier losses on the 
opponent. The more unstable a system gets, the more signif- 
icant will be the effects of initial perturbation which are 
manifested by initial victory and element of surprise. Since 
some systems with large aimed-fire coefficients tend to be 
highly unstable, we can expect this effect to be most 
pronounced in battles involving high-technology and highly- 
lethal weapons. 

•k 

The payoff matrix and optimal probabilities p. and 

k 

g are shown in Table III. The results suggest that the 
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TABLE II 

Effect of Different Perturbations on the Payoff 



Location in 
Figure 5.9 


Co-ordinates of 
Perturbed Point 


Payoff to X 


A 


(6.7, 3.05) 


2.39 


B 


(6.7, 3.025) 


2.41 


C 


(6.7, 2.975) 


2.45 


D 


(6.7, 2.95) 


2.47 




Figure 5.9 Trajectory for Korean War. 



North Koreans should use large initial deployment. In the 
actual war. North KoLea actually deployed almost all of its 
regular force and within a few days captured Seoul, the 
capital cf South Korea. The payoff matrix also shows that 
no matter which strategy is chosen by South Korea, it is 
bound to suffer much more losses than North Korea. Again, 
this is in agreement with history since it is an accepted 
fact that without US intervention, there would be no South 
Korea today. 

So far in the example, we have always considered the 
situation in which the equilibrium point (X,Y) determines 
the replenishment rates as given in equation 5.1. It is 
interesting to investigate the effect on the payoff when the 
initial operating point is at some other location ether than 
<X,Y) . let the new initial point be at (X 1 ,Y 1 ) and consider 
a case where X ! is hept equal to X and only Y 1 is varied. 
(X,Y) has been chosen to be (6.7,3. 0). In Figure 5.10, three 
trajectories corresponding to Y at 2.5, 4.0, 5. 0 are shown 

together with the hyperbolic intersection during stage one. 
Basically, the trajectories correspond to the three stages 
of simulation as before and x is still the victor. However, 
both the payoff (Ly-Lx) and finish time are slightly 
different from operating at (X,Y). Table IV sWws tint y 
inflicts mere losses on x when operating at Y 1 above the 
boundary curve rather than at (X,Y) , but in doing so y is 
defeated faster. Thus depending on his mission, a commander 
can choose to lengthen the battle or inflict more casualties 
on his opponent by choosing a suitable operating point which 
may be other than an equilibrium point. 
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TABLE III 

Payoff Matrix and Mixed Strategies : Korean War 
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YjM6. 7,2.5) 



Figure 5.10 Initial Operating Points at Non-egmilibrium Points. 



TABLE IV 

Effect of Operating at Non-equilibrium Points 



Y 1 


(Ly-Lx) 


Finish Time 
(FINTIM) 


Remarks 


2.5 


2.475 


0.323 


Below boundary 


3.0 


2.435 


0.313 


At equilibrium 


4.0 


2.353 


0.293 


Above boundary 


5.0 


2 o 305 


0.283 


Above boundary 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCIU SION S 

This thesis has covere 
based on the generalized La 
the results has to do with 
the N *N system. The Contin 
be suitable for this pu 
Continuation Methods over n 
and important to our unders 
eguations. The method fin 
accurately and does not 
example to compute the egui 
is presented along with a w 
The derivations and i 
between stability and syste 
portion of the thesis, 
problem, a few interestin 
namely 

(1) Eoth the system 
are intrinsic to a sys 
cf the equilibrium poi 
completely characteriz 
are the differences in 

(2) The dynamics of 
phase trajectories wh 
progresses. Besides 
tion also influence 
curves which separate 
by graphical or backwa 



d a number of subject s which are 
nchester Model. The first part of 
finding the equilibrium points in 
uation Methods have been found to 
rpose. The advantages of the 
umerical techniques are numerous 
tanding of the non-linear set of 
ds all the equilibrium solutions 
need good initial guesses. An 
librium solutions of a 2*2 system 
ay to treat singularity problem, 
nterpretations of the relations 
m parameters form the next major 
By considering the simpler 1*1 
g conclusions have been reached, 

asymptotes and equilibrium points 
tem in equilibrium. The locations 
nts on the x-y and e x , e y planes 

e their stabilities; e and e 

x y 

the system asymptotes ; 
a system are characterized by the 
ich represent the ways a battle 
stability, the domains of attrac- 
the trajectories. The boundary 
these domains can be ascertained 
rd integration. 
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The last portion of the thesis integrates tha concept of 
equilibrium stability and system dynamics. It relates these 
theoretical concept to operational problems. Two operational 
issues are addressed namely , (1) the effect of varying X and 
Y, the initial force deployment on battle outcome, (2) the 
exploitation of stability to influence battle outcome. A 
methodology which combines multistage battle simulation with 
two-person game has been employed and the conclusions are 

(1) Initial force deployment, X and Y can be optimized 
by finding a set of mixed strategies. A suitable pure 
strategy can then be selected from the mixed strat- 
egies ; 

(2) Instability can and should be used to shape the 
course of battle and its outcome. This is particularly 
true in highly unstable warfare which is normally asso- 
ciated with large aimed-fire attritions. Unless defense 
resources are extremely large, it is not possible to 
completely reverse the outcome of a lopsided-battle 
where one side is much stronger than the other. 

As far as military commanders are concerned, the above 
conclusions suggest two things. Firstly, depending on the 
relative strengths, it is not necessarily true that 
deploying the largest possible force will bring victory, 
reduce loss or even buy time. There is an optimum way of 
deploying available forces. Secondly, if a war involves 
large aimed-fire attritions due to weapons like aircraft, 
missiles, tanks, artillery, naval bombardment, etc., then 
initial victory which could perhaps be achieved through a 
preemptive strike certainly affects battle outcome. 
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B. BEC0BZ1EHDATI0HS 



1 . Transf or ma tion of N *N Problem in to the V^_}_ Problem 

It has been mentioned that the simplicity of the 1*1 
problem can be attributed to the simpler mathematics 
involved and our ability to draw and visualize two- 
dimensional pictures. Despite its simplicity, it does share 
many of the properties with the N*N problem. Considering 
these factors, it seems logical that an attempt should be 
made to find the 1*1 equivalent to the N*N system. Another 
reason is that there is much to be gained in terms of 
savings in computational effort by going to the 1*1 
equivalent . 

Of course, the "equivalent" system will not be 
expected to be identical to the N*N problem in every aspect. 
One can only hope that it is equivalent in some sense, for 
example 

(a) Preservation of stability char act er istics and 
dynamics ; 

(b) Preservation of mixed strategies. 

One way of transforming the N*N system into the 
equivalent 1*1 system is t o eguate losses in both systems. 
The equivalent system parameters are obtained by using rela- 
tions such as 



Y 

q e 



EE a ij x iyj 

i j 



where 




coordinates of the equilibrium point in 
the N*N problem 
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eq 



= 



eq 



“ S>: 



3 



The results of the preliminary studies suggest that 
this method of transformation can preserve some stability 
characteristics. The possibility of using the equivalent 
system to obtain the mixed strategies should not be 
dismissed until further studies have been conducted. 

2. Time V aria ble Replenishment Coefficients 

The replenishment rates used in the thesis have 
always been assumed to be constant. In actual wars, constant 
replenishment rates may not be used by either side; at 
times, it may not even be possible to do so. It would there- 
fore be interesting to study the cases which involve time- 
varying replenishment rates r(t). The choice of r (t) , for 
example periodic, non-periodic, ramp, etc. will depend on 
how well it represents practical replenishment rates. 
Whether a mathematical tool can be found to cope with the 
added complexity also has to be considered. 
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C APPENCIX A 

C SOLVING FOP OPTIMUM VECTORS AND 

C 
C 
C 

£********* **** * ******************* ******* ****** ********* 

C THIS PROGRAM COMPUTES THE OPTIMUM FORCE 

C D ISTR1EUTICN, X* GIVEN X,Y ,ANO THE ATTRITION 

C CCEFF IC IENTS. PARAMETERS K1,K2 AND CONFLICT 

C M ATR I > ARE ALSO OBTAINED. 

C THIS IS AN INTERACTIVE PROGRAM. EEF OFE EXECUTION 

C , CHECK TFE VALUES OF ATTRITION COEFFICIENT 

C MATRICES AA,BB,CC,D AND < U 1 , U 2 , V 1 , V2 ) 

C IN THE FILE L PD FORTRAN. THEN ENSURE YOU HAVE 

C ANOTHER FILE CALLEC CCLIB EXEC. TO EXECUTE THE 

C PROGRAM, ENTER "COLIB LPD" AND FOLLOW INSTRUCTIONS . 

C ON INPUT OF X AND Y WILL CREATE THE A MATRIX 

C AND RUN THE LINEAR PROGRAM ZX4LP. 

C***** ******** * <*** *** ************************ ********** 

C 

C A(Ml+M2+2,N*Ml+2) 

C E ( Ml +M2 i 

C C( N) 

C PSOL ( Ni 

C DSOL (M1+M2J 

C RW(I A* l M1«M2*4) + 2*(N*M1) i 

C IW(2*(N + M1J +3) 

REAL A( 4 , 1 J , B (2 ) ,Ci3 ),S,FSUL(3) , DSOL ( 2 J , R W( 34 1 , 

CIER, X,Y,M,PHI ,BB (2 ,3 ) ,D(2,3J 

& ,X1 , X2»Y1,Y2»Y3# K1,K2,CC(2,3J,AA(2,3J ,CO N (5 , 5 J 
INTEGER C,IW< 13) ,M1,M2,IA,N,K 

Q***** ******** ******** ** ***** ********* ************* ***** 

C VARIAEIE DEFINITIONS 

£********* **** ************ ****************************** 

c 

AA.BB ,CC,D ,U1 ,U2, VI, V2=AT TRT I ON COEFF I 01 ENTS 
PSCL, 0SCL=PRI MARY AND DUAL SOLUTION 
M1=NUMEER OF EQUALITY CONSTRAINTS 
M2=N UME ER OF INEQUALITY CONSTRAINTS 
RW,IW=PARAMETERS USED BY IMSL RUUTINE ZX4LP 
I EB= ERROR CODE FROM ZX4LP 
N = NUM E E R OF UNKNOWNS 

A=MAT R I X DEFINED IN EQUATION A. 3 OF PAPER 
"LANCHESTER EQUATIONS AND GAME THEORY" 

BY P.H.MOJSE ANCJ. M. WOZENCRAFT 
C- VECT CR CQNT A I NI NG COEEFICIENTS OF 



C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

E 

C 

C 

c 



UF OBJECTIVE FUNCTION IN LP „ _ 

B=VECi OR CONTAINING THE RIGHTHAND SIDES 
OF THE CONSTRAINTS. 

I A=M1+M2+2;IS THE ROW DIMENSION OF A MATRIX 
Q = CQL UMN DIMENSION OF A 



Ml =2 




M2=0 




N=3 




IA=M 1 + M2 + 2 


C = N+ Ml+2 


P=M1 +M2 


Ui = 


.30 


U2 = 


.30 


VI = 


.10 


V2 = 


.0345 


V3 = 


2.0 


AA (1 


,1) = 1.0 


A A ( 1 


,2) = .3 


AA ( 1 


,3) = 0.0 
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A A ( 2 , 1) = 


.60 


A A (2 , 2 ) = 


.900 


AA (2 ,3) = 


0.0 


EE(ltl) = 


0.15 


BB (1 ,2) = 


.10 


E B ( 1 ,3) = 


0.0 


£B (2 , 1 ) = 


0.15 


BB ( 2 ,2) = 


.3 


EB (2 ,3) = 


0.0 


CC(ltl) = 


1.2 


CC(1, 2) = 


1.0 


CC(1, 3) = 


15. 


ec ( 2 , i ) = 


1.1 


CC (2 ,2) = 


0.6 


CC (2 ,3) = 


15.0 


0(1,1) = 


.00 


0( 1 , 2 J = 


.00 


0(1,31 = 


0.00 


0(2,1) = 


.00 


0(2,2) = 


.00 


0(2,3) = 


0.0 



£***#* **«*<**« * * **$$:(:£***$ *$***$$ ***** ********■<■ ********* 
C REPET I7ICN LOOP STARTS HERE 

£***** ***# Jf $ $*$*** **>£ * **•$** **# S********* 

c 

CCNT INUE 



7000 

C 

C 

C 

296 

C 

C 



ENTER X AND Y 
W R I T E ( 6 » 29 6) 

FORMAT (IX, 'ENTER X ANO Y 
R EAD ( 5 ,* ) X, Y 



B 1= B E ( 1 , 1 J+B B ( 2, 1) 

62= BE (1 ,2) + 6B( 2, 2) 

B 3= BE (1,3)+8B(2, 3) 

01= D ( 1 » 1 ) +0 ( 1,2) + D C 1,3} 
02= 0 (2 ,1)+D( 2,2) +C( 2 , 3) 



ONE AT A TIME') 



£ 

& 

& 

£ 

& 



999 

C 



1000 

C 



A (1, 2)= CC<1, 2) 
/ XU5 . 

A (1*3)= CC(1 , 3) 
/ X+l 5 • 



A (2, 1 ) 

/ X + l 5 • 



CC (2, 1) - 



A 12.2)= CC (2,21 - 
/ X + l 5 « 

A (2,3)= CC (2, 3) - 



B ( 1) = 1 .0 
B (2) = 1.0 
CONTINUE 

C (1) = 1.0 
C (21 = 1.0 
C (31=1.0 



AA ( 1 , 1 ) 


+ 


( Dl-Ull/Y 


• 


( Bl-Vl) 


AA ( 1 ,2 ) 


+ 


( ci-un/Y 


- 


( B2-V2 ) 


AA ( 1 ,3 ) 


+ 


( Dl-UD/Y 


- 


( B3-V3) 


AA ( 2 , 1 ) 


♦ 


( C2-U2 )/Y 


- 


(Bl-Vl) 


AA ( 2 ,2) 


+ 


( D2-U21/Y 


- 


( B2-V2) 


AA (2 ,3 1 


+ 


( D2-U2 )/Y 


- 


( B3-V3) 


, A ( I » J ) 
,13,13, 




1 ,2) ,1=1, 
' »F10.3) 


2 ) 
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130 



C 

200 

300 

c 



WRITE ( ] 4 » 130 J 

FORMAT (IX, 'COEFFICIENTS OF CONSTRAINTS (AJ») 



DO 30 C 1 = 
WRITE ( 
TH 

FORMAT 
C QNT i NUE 



»R 

A, 2 00 J ( A ( I ,KJ ,K=1 »NJ 

S NUMBER = N+Ml+2 
IX, 7F 15-6 J 



CALL Z X ALP (A , I A , B ,C, N ,M1 , M2, S ,PSUL ,DSQl,RW , 
I Wil EF i 



£************** 4$********* ****** ****** ****************** 

C CCMPLTE VI , Y2 ,Y3 ,X1, X2, H 

C***** **** *$£* ******** **** * ****** ***** ****** ** ********** 

c 

P H I = S 

Y l = PSCl( 1J*Y/PHI 

Y 2 = PS C L l 2 ) *Y / PHI 

Y3=PSCL (3) *Y/ PHI 
X 1=DS L L ( 1 J *X/ PHI 
X2-DSCL (2)*X/ PHI 
M = X* Y / (PHI-15 1 
K 1= ( <E1-V1I*Y1 4 

K2=( (C1-U1J*X1 4 



(B2-V2MY2 4 (B3-V3J*Y3 •* M J/X 
(D2-U2 ) *X2 - M J/Y 



C * ** * * * * ** *** ************* *********** ******************* 
C COMPUTE CONFLICT MATRIX 

C* ************* 4 ** * *** ********* ************************* 

c 

C ON ( 1 , 1 J = LI + AA ( 1 , 1 ] *Y 1 + AA(1,2J*Y2 + AA(1,3J 
6 *Y3 

C ON ( 1 , 2 i = 0.0 

C ON ( 1 ,21= AA(1,1J*X1 4 BB ( 1 , 1 J 
C0N(1,AJ = AA(1,2)*X1 ■* BB ( 1 , 2 J 
C ON ( 1 , 2 J = AA ( 1 , 3 J * X 1 ♦ BB < 1 , 3 ) 

C 0N( 2 » 1 )= C.O 

CON(2,2J = U2 ♦ AA ( 2 » 1 ) *Y 1 + AA(2,2J*Y2 + AA(2,3J 
C * Y3 

CON(2,2) = AA(2,1)*X2 4 BB(2,1J 
C QN( 2 , A J = AA(2,2J*X2 4 BB(2,2J 
C 0N( 2 , 2 J = AA ( 2,3) *X2 4 BB(2,3J 
C ON ( 3 , 1 i= CC( 1 , 1 J ♦ Y1 4 0(1. 1J 
C 0N( 3 , 2 ) = CC(2,1J*Y1 4 0(2. 1J 
C ON ( 3 , 2 ) = VI 4 CC(l,li*Xl 4 CC(2,1)*X2 
C ON ( 3 , A I - C.O 
C QN( 3 , 5 i = C.O 

C CN ( 4 , 1 ) — CC(1,2J*Y2 4 0( 1, 2J 
C ON ( 4 , 2 J - CC(2,2)*Y2 4 D(2,2) 

C 0N( 4 , 2 )= C.O 

C ON ( 4 , A J = V2 4 CC(1,21*X1 4 CC(2,2J*X2 
C 0N( 4 , 5 I - C.O 

C ON ( 5 , 1 J = CC(1,3J*Y3 4 D(l,3) 

C 0N( 5 , 2 ) = CC ( 2 ,3) *Y3 4 0(2, 3J 
C 0N( 5 ,21= 0.0 
C 0N( 5 , A i = C.O 

C 0N( 5 ,5 i= V3 4 CC ( 1 , 3 J *X 1 4 CC(2,3J*X2 
DO 2 1 C C LA = 1 , 5 
CO 2110 L5= 1,5 

C^N ( LA,L 5 J = -1.0*CON(L4, L5J 



CONT 
CGNT INUE 



2110 
2100 
C 

£************* ****************************************** 
C PRINT RESULTS 

q* ****** ****** ****************************************** 
C 

WRITE(lA,150i 
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150 FORMAT (IX. 'RIGHT HANO SIDES OF CONSTRAINTS ( B i 4 ) 

WRITE ( 14, 160) (B(I),I=1,R) 

C THIS NUMBER =M1+M2 

160 FORMAT (IX ,2F1 0.4) 

WRITE ( 14 ,170) 

WRITE < 14 , 180) ( C( 1 ) , I = 1, N ) 

170 FORMAT (IX, 'COEFFI CIENTS OF OBJECTIVE FUNCTION(C)' 

C THIS NUMBER =N 

180 FORMAT <1X,3F10.4) 

WRIT E ( 14 , ISO! N 
WRITE (14 ,1511 Ml 
WRITE (14 ,1521 M2 

190 FORMAT ( IX , 'NUMBER CF UNKNOWNS (N1 = ',1201 

191 FORMAT (IX, 'NUMBER CF INEQUALITY CONST RAINT S( Ml)=' 

& ,151 

192 FORMAT ( IX , 'NUMBER OF EQUALITY CGNSTRA I NTS ( M2) = • 

£ ,171 

WRITE ( 14 ,210) S 

210 FORMAT (IX, 'VALUE CF OBJECTIVE FUNCTION(S) =' 

£ , F 15 . 6 1 

WRITE ( 14 ,220) ( PSCL( I) f I=l,N) 

C THIS NUMBER = N 

220 FORMAT (IX, 'VALUE CF PRIMAL SOLUTION (PSOL) = ' 

£ , 3F1 5 • 6 ) 

WRITE (14 ,230) ( DSCLi I ) , 1=1 ,R ) 

C THIS NUMBER = M1+M2 

230 FORMAT (IX, 'VALUE CF OUAL SOLUTION (DSOL) = ' 

£ , 2F1 5 . 6 ) 

WRITE ( 14 ,240) I ER 
240 FORMAT (IX, 'IER = • ,151 

WRITE ( 14 , 570) X 
WRITE 04,580) Y 
WRITE (14 ,520) XI 
WRITE ( 14 ,540) X2 
WRITE ( 14 ,500) Y 1 
WRITE ( 14 ,510) Y 2 
WRITE ( 14, 520) Y3 
WRITE ( 14 ,550) M 
WRITE ( 14 ,660) K1 
WRITE ( 14,670) K2 



570 


FORMAT (IX, 


•X 


— 


fl 


,F15.6) 


580 


FORMAT (IX, 


•Y 


SL 


1 


,F1 5*6) 


500 


FORMAT (IX, 


'Y 1 


= 


fl 


,F15.6) 


510 


F CRM A T ( IX , 


•Y 2 


S, 


fl 


, Fl 5 • 6) 


520 


FORMAT (IX, 


•Y 3 




fl 


, F 1 5. 6) 


530 


FORMAT (IX, 


•XI 


•= 


fl 


,F15.6) 


540 


FORMAT (IX, 


•X2 




fl 


,F15.6) 


550 


FORMAT (IX , 


•M 


=: 


fl 


, Fl 5 *6) 


660 


FORMAT (IX, 


•K 1 




fl 


, Fl 5 .6) 


670 


FORMAT ( IX, 


•K 2 


= 


fl 


,F15.6) 



C 

C***** ******* ********* **************** ******** ********** 

C PRINT PARAMETERS 

Q* ****** ****** ******** ********************* ************* 

c 



WRITE ( 14 ,705) 
WRITE ( 14 ,710) 
WRIT E ( 14, 715) 
WP1TE(14,7 20 ) 
WRITE ( 14,725) 
WRITE ( 14 , 7 £9 ) 
WRITE (14 ,750) 
WRITE ( 14,7 70) 
WRITE ( 14 ,750) 
WRITE (14,7 80) 
WRITE ( 14,759) 
WRITE (14,750) 



U1 

U2 

VI 

V2 

V3 



( ( A A ( I ,J) ,J = 1,3 ) ,1 = 1,2) 
( (BE ( I , J ) , J=l, 3 ) ,1=1, 2) 
( ( C C ( I , JJ ,J = 1,3 ) ,1 = 1 ,2) 



) 
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WFITE (14,7 


50) ((0(1,J),J-=1,3) 


,1=1 


, 2 


70 5 


FORMA! (IX, 


'U1 = ' 


»F1 5 .6) 






710 


FORMAT ( IX , 


•U2 = • 


,F 1 5 . 6) 






715 


FORMAT (IX, 


'VI = • 


,F15 .6) 






72 0 


FORMAT (IX, 


' V2 = ' 


,F1 5.6) 






725 


FORMAT (IX, 


' V3 = ' 


,F1 5 .6) 






750 


FORMAT (IX , 


2F1 5.5) 








760 


FORMAT ( IX, 


'VALUE 


OF A MATRIX 


(AA) 


' ) 


770 


F ORMAT ( IX, 


'VALUE 


OF B MATRIX 


(BB) 


' ) 


780 


FORMAT (IX, 


•VALUE 


OF C MATRIX 


(CO 


• ) 


790 


FORMAT (IX, 


'VA LUE 


OF 0 MATRIX 


(0) ' 


) 



c 

C*****XXXXX*XXXXXXXXXXXXXXX XXX* XX XXXXX XXX XXX* XXXXX* X**** 

C PRINT CCNPLICT MATRIX 

C***XX XXX* XXXX XXX*XXX******X**X**X***XX****XX*XX**XX*XX* 

c 



WRITE < 14 ,673) 

WRITE (14,674) ( <CCN< L4.L5) ,L 5= 1,5) .14=1,5 ) 

673 FORMAT (IX, ‘VALUE OF Th£ CONFLICT MATRIX (CON)*) 

674 FORMAT (lX,5F14.o) 

C 

C**XXX XXXXXXXXX XXXXXXXXXXX XXX XX XX XXXXXXXXXXX XX XXXXXXXXXX 

C PRINT MEANING OF I ER 

C*XXX*X*XXXXXX XXXXXXXX XX XX ***4 XXX XX XXX XXX X XXXXXXXXX XXXXX 

C 



390 

400 

410 

420 

430 

440 

450 

460 
C 



£ 

& 

£ 

£ 

£ 



WRITE ( 14,350) 

WRITE ( 14 »4C0 ) 

WRITE ( 14 ,410) 

WRITE < 14 ,420) 

WRITE ( 14 ,420) 

WRITE ( 14,440) 

WRIT E ( 14 ,450) 

WRITE 114,460) 
FORMAT 1 IX, • • ) 
FORMAT (IX ,*IER-130 
FORMAT 1 IX, 'IER=13 1 
• ) 

FORMAT (IX, 'IER= 132 
CONSTRAINTS' ) 
FORMAT (IX , ' I ER= 13 4 
UNBOLNCED* ) 

FORMAT (IX, »IER=135 
I NFEASIELE') 

F ORMAT (IX. 'I ER= 13 6 
SOLUTIONS* ) 

F CRM AT (9X , ' DO NOT 



INCICATES 

INCICATES 

INCICATES 

INCICATES 

INDICA TES 

INCICATES 

SATISFY THE 



M2>N ' ) 

EXCESSIVE ITERATIONS 
REOUNCAN C IES IN 
OBJECTIVE FUNCTION 
CCNSTRAI NTS 
PRIMARY GR DUAL 
CONSTRAINTS' ) 



C****X****XX**XXXXXXX* X s X XX xxx x xxx* xxxx xxxxxxxxxxxxxxxxxx 
C REITERATION ROUTINE 

C**XX* XXXXXXXX X X XXX XX XXX XX XXX XX XXX XXX X XXXXXXXXXXXXXXXXXX 

c 



WRITE (6,470) 

WRITE (6,480) 

470 FORMAT (IX, 'TO CONTINUE ENTER 1') 

480 FORMAT ( IX, 'TO STOP ENTER -1') 

R E AD ( 5 » * ) L7 

lF(ABS(L7*1.0).LT.1.0E-5) GOTO 9000 
WRIT E ( 14 ,450) 

WRITE (14,491) 

WRITE ( 14,452) 

490 FORMAT (IX,' ' ) 

49 1 F QRM AT (IX,'** **** x**xx**xx***xx*x * xxx x x**xx* • ) 

492 FORMAT (IX, • • ) 

GOTO 70 0 C 

9000 CONTINUE 
STOP 
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APPENDIX B 

NUMBER OF EQUILIBRIUM SOLUTIONS OF THE N*N PROBLEM 

W € first consider the 2*2 and 3*3 problem and extend the 
results to the N *N problem. 

1 • 2*2 Problem 

In general / each trivial solution leads to a unique 
equilibrium solution in a continuation process and there 
will be egual number of equilibrium and trivial solutions. 
Section 4 discusses what happens when there is degeneracy. 
The trivial solutions are obtained by solving 

Z l (U l +a ll 2 3 +a 12 Z 4 5 = 0 
z 2 (u 2 +a 21 Z 3 +a 22 Z 4 ) = 0 
z 3 (u 3 + C ll Z l' fC 21 Z 2^ ) = 0 
Z 4 (U 4 + C 12 Z 1 + C 22 Z 2 } . = 0 

At first glance, there would 

solutions corresponding to the number 

lefthand sides of equation B. 1 zero. Each lefthand sides can 

be made zero by either making z^ or the terms in parenthesis 

equal to zero. But closer examination reveals only six 

allowed cases, in non-degenerate cases, corresponding to (1 
2 

+2 + 1) =6 solutions. Table V shows how these cases 

arise. 



(eqn B.l) 



seem to be 2 4 trivial 
of ways of making the 
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TABLE V 

Trivial Solution for 2*2 Probem 



Case z. which are made zero 
1 



Remarks 



1 

2 

3 

4 

5 

6 
7 





4 cases 



all terms in parenthesis 
equal to 0 
degenerate case; 
see section 4 . 



2 « 3*3 pr ob lem 

Here, the trivial solutions are obtained from 

z l (u l +a ll Z 4 +a 12 z S +a 13 z 6> = 0 
z 2 ( u 2 + a 2i z 4 + a 22 z 5 + a 23 Z 6- > * 0 
z 3^ u 3 + a 31 Z 4 + a 32 z 5 + a 33 z 6'* = 0 
Z 4^ U 4 +C 11 Z 1 +C 21 Z 2 +C 31 Z 3^ = 0 
Z 5 (U 5 + C 12 Z 1 + C 22 Z 2 + C 32 Z 3 :I = 0 
Z 6 (U 6 +C 13 Z 1 +C 23 Z 2 +C 33 Z 3 ) * 0 

Again, there are some cases which are not allowed 
because cf inconsistencies in non-degenerate cases. Table VI 
shows the different cases. 

In this case, the total number of allowed cases is 
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1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 



TABLE VI 



Trivial Solution for 3-3 Problem 



z. v;hich are made zero Remarks 

l 









9 cases 



all items in parenthesis 
equals to 0 . 

degenerate case; 
see section 4. 

degenerate case; 
see section 4. 
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3 . Exte nsio n to N* N Problem 

Now we know the type of allowed cases, it is easy to 
recognise the pattern of results. The pattern looks like a 
Pascal triangle with each element squared. 



N*N 


, 2 


2 


1*1 


1 + 


1 


2*2 


l 2 ♦ 2 2 


♦ l 2 


3*3 


l 2 + 3 2 + 


3 2 + l 2 


4*4 


l 2 + 4 2 + 6 2 


+ 4 2 + l 2 


5*5 


l 2 + S 2 + 10 2 + 


10 2 + 5 2 + 



Number of 
Equilibria (NjO 

= 2 
= 6 
= 20 
= 70 
= 252 



For N*N in general. 





2 



i=0 



N can also be written as 




The proof of the identity. 




i = 0 

can be found in [Ref. 12]. 
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4 . Deg ene racies 



For the 2*2 problem , a degeneracy can arise when z 1 
(or any ether z ) is made equal to zero. The other three 
lefthand sides in equation B. 1 are made equal to zero making 
the terms in the parentheses egual to zero. 

In general, such a case does not correspond to a new 
trivial solution because there is an inconsistency when z 0 = 
-u 3 /c and z 2 = -u 4 /c 22 . However, it seems that by making 
u 3/ c 2 1 = u 4 /c 22 * inconsistency no longer exists and 

there will be an infinite number of trivial solution as long 
as z 3 , z 4 are chosen to satisfy u 2 + a 21 z + a 2 2 Z 4 = ® * it 
turns out that the number of equilibrium solutions still 
remains at a maximum of six (disregarding the case with 
infinite number of equilibria). Similarly for the 3*3 case, 
even if there is degeneracy, the number of equilibria will 
not exceed 20 . 

That the above is true can be shown by considering 
the 1*1 problem. In this case, the two non-degenerate cases 
correspond to (a) x = y =0 (b) v + cx = 0 and u + ay = 0. A 
degeneracy can occur if v = 0 or u = 0 in which case there 
seems to be an infinite number of trivial solutions lying on 
the y or x axis respectively and hence an infinite number of 
equilibria. But when the actual hyperbolas are plotted, 
there are only two intersections and hence two equilibrium 
points. Furthermore, when the Continuation method is used to 
find the equilibria, there are only two equilibria irrespec- 
tive of whether the trivial solutions are chosen to be 
degenerate or not. 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c****« 

c 



APPENDIX C 

P RGGR A M LISTING F UR SOLVING 2*2 
USING CONTINUATION METHOD 



SYSTEM 



£ 

£ 

C 

C 

c***** 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c****« 

c 

c 

c 

c 



THIS 

ECU! 



REAL 
XT ,D 
PEAL 
INTE 
#CGN 
, FAS 
DIME 
XA (1 
,XD( 
EXTE 
COMM 
COMM 
COMM 
fOMA 
DATA 
MM/1 
ERFL 



*** * 


* **** *** 


**** 


PRO 


GFAM SOL 


VES 


LIBR 


IUM PCIN 


TS U 


ICN 


1 11 8 FUR 


THE 


**** 


* *$:**.$$:£ 


**** 


* 8 


A . X , CE L 


t T , W 


X ,QMAX , XFTEM 


r , TL 


* 4 


XA ,XE,X 


C , XD 


GER 


I i J , K ,N , 


M , MM 


, KK , 


NRECGN, S 


IGN 


T , QC 


OUNT .RUN 


. ERF 


NSI 0 


N WKAREA 


i 100 


0C5 i 


,XE (1C05 


J ,XC 


ICO 5 


i ,TA( 100 


5) ,D 


RNAL 


f CN 




ON / 


PEPLEN/ 


R 14) 


ON / 


PARAM/ A 


(2,2 


ON / 


TREPAR/ 


XT (4 


X »QC 


OUNT ,FUN 


, XFT 


DEL 


/ l.OD-16 


/ » 


/ ,M/ 


1/ , NS 1G/ 


5/, I 


AG/O 


/, REFLAG/O/ 


**** 


******** 


**** 


A EL E 


CEF INAT 


IONS 


**** 


* * * *c* *** 


**** 


EE E 


GUATION 


3-7 


,S AB 


CVt 




E RGCT OF TH 


E 2* 


ONTI 


NUAT ION 


PROC 


RAME 


TER DEFI 


NING 


1C 1 . 


C 




SMAL 


L TIME I 


NT ER 


PART 


1AL TIME 


DER 


$,C$ 


»DS = T HE 


0R1G 




LANC 


HEST 


CCRR 


ESPOND T 


0 TH 


MENT 


CCEEF. 


IN T 


ART I 


T JON INT 


ERVA 



A 2*2 SYSTEM FOR THEIR 

SING CONTINUATION METHODS. SEE 

CRY AND IMPLEMENTATIONS. 

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * ** 

KAREA.FNORM, WK,U,R,B,C,D,EP,AP, 
EFT, TS,F1,F2 ,F3,F4,ERP 
» TA 

»COUNT»IA,IDGT ,IER,NSIG» I E, I REP 

LAG,REFLAGfCHFLAG»CTR ,CP1 ,JJ 
J ,WK( 11 5 J » BP ( 4 J , APA4.AJ , 

41005 J 
X ( 4 J 



) , B ( 2 ,2 i , C ( 2 , 2 J » D (2 , 2 ) » U ( 4 J 
) ,X (4 j , T (4J , TS, FAST, CON » COUNT ,N 
EMP , SUM, SIGN , IER 

T MAX/200 / , IDGT/O / ,IA/4/ , 



THE MEANING. 

STEM TO BE EVALUATED EY THE 
HOMGTQP Y JHAS VALUE BETWEEN 
APPROXIMATE THE 
PARAMETERS IN THE 



TO 



CONTINUATION PROCESS 
COUNT^COUNTER FOR TH 
THE CGNTINUATI 
CON= THE CCNDITION FO 
FCN= A SUBROUTINE USE 
WK,NSIG , ITMAX ,IER,FN 



VAL USED 
I VATIVE. 

INAL A,B,C,0 
ER ECUATIUN 
E SELF ATTRI 
HE LANCHESTE 
l FRCM T -0 T 



TION £ REFLENISH- 
R EQUATION. 

C T = 1 IN THE 



IE ,1 A=PARAMETERS IN 
CTR= CQUN T ER FOR PLOT 
NRECON^REPE TITION C 



E # OF TIME STEPS ADVANCED IN 
ON PROCESS. 

R REACHING T 
C BY THE IMS 
GRM=PARAMETE 
ZSCNT. 

THE ROUTINE 
TING THE CUR 
GUNTER. 



= 1 

L ROUTINE ZSCNT. 
RS IN THE ROUTINE 

LEQT IF. 

VES 



* ******4 ************ ************ ***** * ** * ********** 

INPUT ATTRITION CCFFICIENTS, COUNTERS AND FLAGS. 

**** **** 4 ** ****** *********** * **** ************** 



All, 1)=1 .000 
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a( i 

A(2 
A( 2 
E ( 1 
E(1 
E( 2 
E( 2 
C(1 
C(1 
C ( 2 
Cl 2 
D( 1 
0(1 
0(2 



2 1= C . 3D0 
11 = 0 .600 
2 i = 0 .900 
11 = 0. 1500 
21=0.300 
11=0. 1500 
2 1=0.300 
11 = 1 .2000 
21 = 1 .000 
11 = 1 .1000 
21 = 0. ECO 
11=0 .OCDO 
21 = 0.000 
11 = 0.0000 



0(2. 21 = 0. 00 C 
U(l) =0.300 
U( 21 =0. 3 DO 
0(31 =0.100 
Ul 41 =C. 200 
N=4 

T( 11 =0.000 
1(21 =C. 000 
T ( 31 =C. 0 DO 
1(41 =0.000 
C C UN 1 = 1 
CC0UNT=0 
PGN= 1 
I ER= 0 
1R EP = 1 
NREC 0N=0 
C 

C CALCULATE P VECTOR USING ONE POSITIVE EQUILIBRIUM 

C SOLUTION OR MODIFY R VECTOR AS APPROPRIATE. 

*$**•#*** * 4 ********************* ********* **** ** 

c 

X ( 11 =0. 6 1 5 3 6 DO 
X(21 =0. 384600 
X( 31 =3.0 7 69DC 
* (41 =0.923100 

R(1)=X(1J 3 MU(1J + A(1,11*X(31+A(1,21*X(411+B(1,11 
G *X (3 1+B( 1 ,21*X(41 

F<21 =X( 21*(U(21+A(2, 11*X(31+A<2, 21*X(411 +B(2, 11 
& *X (3 1+B( 2 , 2 1 * X( 4 1 

P(31=X(3 1*(U(3J+C(1, 11*X(11+C(2, U*X(211+Cll,ll 
G*X(1 1+D( 2 ,11*X(2 1 

P(41=X(4)*(U(41+C(1,21*X(11+C(2,21*X(2J1+0(1,21 
£*X(1 )+D(2 ,2 1*X(2 J 
WRITEU .9998) (I, X(I1 ,I.R( II , 1=1,41 
9998 FQRM/iT(*0*f3X, , X( , ,l3,'J = , ,rl2.4,4X, ' R ( 1 » 13 » ' 1 = ' . 

Of 12.41 
C 

c 

c WHILE IREF (TFE CODE FOR REPEATINT 1 IS 1 .THE CGNT- 

C INUATICN FRCCESS WILL BE REPEAT EC-EACH REPE TITIGN 

C WILL INCREASE TS IN MULTIPLES OF 0.005. 

C 

**** ^fr#*****##^# ***#£*$$# $ $*$$$#:$ 

C 

7 IF(. NOT. (IREP.EQ .111 GO TO 8 

NRECCN=NPEC0N+1 
C 

c***** ** *•£££$:{»$ **£*$****$#:£ 

C SET TRIVIAL SOLUTION 

£#**$$**** ************ ******************* 

c 
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X ( 1) = -C. 10C00D0/1 . 200D0 
X (2J=-C*0003D0/1. 100D0 
X (3) =-C.30C0/ 1.00000 
X (4)=-C.C000/0. 30000 

t m=c.cco 

T ( 2) = C . CDO 
T (3) =0 .OCO 
T ( 4) = C . CDO 
C HFL A G S G 
QCOUNT-O 
R 0N= 1 
ERFLAG=C 
R EFLAC--G 
C 0 UN T = 1 
I ER=0 
C TR= 1 

TS=DFLCATCNRECGN) *0.005000 
C 0N= I D 1NT ( ( 1 . OCO-TaJ )/TS J +2 
C 

£************* =#******: ****** ** ******** ************** ****** 
C VnHIL E T<1 or corrector step ooes not return with 

C ERRUR FLAG IER 

C***** ******** * ************** ****************** ********** 

c 

5 IF(.N07.((IER.NE. 129) .AND.(1ER.NE.130) .ANO. 

0 ( IER. NE. 131) .AND. (COUNT. LE.CONJ) j GO TO 6 

C 

C TREAT THE SINGULAR CASE IF REFLAG OR ERFLAC- 

C IS SET. 

C 

CALL TREAT (REFLAG, ERFLAG,CHFLAG) 

C 

£***** ******** 4 ************** ********* *************4 ***** 

C FREDICTOR ''TEPJ 

C SET UP THE~SET*OF DIFFERENTIAL ECUATIGN ANO 

C CALL LE G I IF TO SOLVE THE SYSTEM AX=B 

q***** ***************************************** ********** 
C 

DO 101 1=1, MM 

EP(1)=( R( 1)-B( 1» l)*X(3)-6( 1 ,2)*X (4) ) *DEL 
EP ( 2 )= ( R(2)-E(2»1)*X(3)-B( 2 » 2 ) *X (4 ) ) *DE L 
BF(3)=( R(3)-Dl 1,1)*X(1) -0(2, 1)*X (2) ) *JEL 
E F <4 )= ( R(4)-0( 1,2)*X(1)-D(2,2)*X(2))*DEL 
AP (1 ,1) — U 11 )+A( I f l)*X (3 )+A( i ,2) *X( 4) 

A F ( 1 ,2 ) =0 .0 CO 

AP(1 ,3)=A(1 f l)*X(l)*8(l ,1)*T(1) 

AP(1,4) =A(1 ,2)*X(1)+B(1 ,21*7(1) 

AP ( 2 ,1 ) =0 .0 CO 

AP (2,2) =U (2 )+A( 2,1) *X( 3 )+A( 2,2) *X( 4) 
AP(2»3)=A(2»l)*X(2)+6(2 ,1)*T(2) 

AP (2 ,4) =A(2 ,2)*X12) + B(2 ,2) *7 (2) 
AP(3,l)=C(ltl)*X(3)*0(l ,1)*T(3) 

AP(3,2)=C(2 ,1)*X(3)+D(2 » 1 J *T (3) 

AP(3,3)=U(3 )+C(l,l)*X(l )+C(2,l)*X(2) 

AP ( 3 ,4) =O.OCO 

4P(4,1)=C(1 t2)*X(4MD(l ,2)*T(4) 

AP(4,2)=C(2 ,2)*X(4) *0(2 ,2)*T(4) 

AP ( 4 ,3) =0.0C0 

AP(4,4)=U(4)*C( 1,2)*X(1 ) + C( 2,2) *X(2) 

CALL LE GT IF ( AP ,M,N, I A, BP, IDuf,W RARE A ,IE) 

CO 102 K= 1, N 

X (K) =BP( K)+X(K) 

102 CONTINUE 

101 CONTINUE 

C 

£************** ******* *********************************** 
C CORRECTOR STEF: 
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4 $ *afca{:$ **# *a££ afrJfrafraicfcJfrafrafraitafrataJtaCcCt afrafcajt £££ 

c 



CALL ZSCNT ( FC N ,NS J G, N , ITM AX , T , X ,F N CRM , WK , I ERJ 
UM A X=2 5C. 0D0*DM IN 1(DX( 1J,DX(2J ,DX(31,DX(4jj 
C 

C COMPUTE DIFFERENCE BETWEEN PREVICUS AND CURRENT X 

C 

DO 104 K=1 , N 

CX(K 1=D ABS( XT( K 1-X ( K 1 1 
104 CONTINUE 

C 

C CALL SUBROUTINE TO DETECT FAST CHANGING COMPONENT 

C 

CALL DETECT (DX , REFLAG* ERF LAG) 

C 

C 

£###$$ ***# * 4 $ jCt#*******:*'#****#^## 

C COMPUTE ERROR IF T IS NEAR 1.0 

C********* :$.$*$**$*:$<*:£*$*# S*#***^# *****=x*^t* af $a{c:Je 

c 



IF ( .NOT • ( DABS( T (1J-1.0D0) .LE. 0.000 5 DO) 1 GO TO 
£ 1031 

f 1=-X( 1 1*(U (1)+A(1, 1)*X (31+A (1, 21*X(41 1 
£ «R<11-B(1,1 J*X(3)-B(1,2 )*X(4) 

F 2=-X( 2 1 * (U (2) + A (2 , ll^X (31 + A (2 , 2 1*X ( 41 J 
£ +R(21-B(2»1 1*X(31-B(2,21*X(4) 

F3=-X( 3 1*(U (3 1 +C (1 » 1 1*X ( 1 1+C (2, 1 1*X ( 21 1 
£ +R(31-Di 1,1 1*X( 11-0(2,1 1*X( 21 

F4=-X(4 )*(U(41+Cll,21*X(ll+C(2,21*X(2Jl 
£ + R(4)-D(1,21*X<11-D(2,21*X(2) 

ERR=CSQ RT (F 1**2+F2**2+F 3**2 + F4**21 
WPITE(1 ,997 1T( 1 j , ERR 
1031 CONTINUE 

C 

C OUTPUT VALUE OF X VECTOR WHEN T IS NEAR 1.0 

C***** $*** * ** *4***#$** #ajt afcajt * afrjfrafc* afr 

c 

IF(.NGT.(T ( 1 1 . G E. (0.99 DO) J I GO TC 863 
WRITEd ,935 1T(11 
fcRITE( 1 ,999 HK,X(K1 ,K=1 ,N) 

865 CONTINUE 

C 

Qa { 14444444444444 4444444 #444 444 4 4 44 44444 *$**#$* 4 * 4 4 44 * 44 

C CHANGE DOUBLE TO SINGLE PRECISION FOR PLOTTING 

Q44444 4444 Jfrajsajta* ^aGa&^afraCufraJtafrJfr** * 4 a* * * **a(iaCtaJi*ajt * *a{c * * * 44 $ * V * * J****** 

c 



6 

C 



TA ( CTR 1 =SN GL (T ( 11 1 
XA (CTR) -SNGL (X ( ll 1 
XB (CTR)=SNGL(X (21 1 
XC (CTR1=SNGL(X (31 1 
XD (CTR J =SNGL (X (41 1 
CQUNT=CCUN T+ 1 
CTF^CTR+l 
GO TO 5 
CONTINUE 



Q44444 4 444 4 444 4 *444444 4 444 444 444444444 

C COMMENCE PLOTTING ROUTINE 

£:((** ajcagc af **afc 444 afc * at * afr * ** * a{t* 44 *■$:{: * a{t*a(sa{i **a{t * ajtajt* afr * a} 44 

C 

CALL TEK618 
C CALL CONPRS 

CALL PAGE (14. 0, 10 .51 
CALL NOBRDP 
CALL B L OWU P( 0.41 
CALL AFEA2D(12.0, 8.01 



4 4444 444444 
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CALL XNAMECTIME IN SECSS',1001 
C ALL Y NAME (* X(I1 , 1 0 0 J 

CALL XNONUM 
CALL YNONUM 
C 

C INPUT HEADING IF REQUIREC 

C 

C CALL HEADIN( ' INPUT HEADING HE RE**** ******* ****** 

C ,1C0,1.,41 

C CALL HEADINC * ***********AND HERE*** ******** ****** 

C ,100,1., A) 

C CALL HEADIN( • *********AND HERE******** • ,100,1. »*tl 

C CALL HEADIN( • *********AND HER E**** *** * *#**<£$ • 

C ,100,1., A) 

CALL GFAF (0.0 , * SC ALE* , 1.00,-56.0, * SCALE* ,35.01 
CP1=CTF-1 

CALL CUKVECTA ,XA, CPI ,-l) 

CALL CURVE (TA,XB, CPI, -11 
CALL CURVE (TA ,XC, CPI ,-11 
CALL CURVEtTA , XO, CPI , 1 ) 

CALL ENOPLIOl 
WRITE U ,99101 

9910 FORMAT (*0* ,3X ,* ENTER 1 TO CONTINUE * ,3X , 'OTHER NO. 
£ TO STOP* 1 

RE AD (5, *1 IREP 
GO TO 7 

8 CCNTINUE 

CALL CONERL 
STOP 

955 FORMAT( ' 0 * , 3X ,' T = ' ,D2A.12) 

99 7 FORM AT( ' 0 • , 2X ,* T - • »D24.14,3X,'EkRGR=* ,02 A. 61 

99 9 FORMAT! • O' » 3X »' X ( • ,13,' 1=*,D24.14/) 

END 

at afra(cafr at ^ at afrafr atajtat afc afloat at =ta{s * at atafc * atafrafc * at at* at at at at at at atatatat*** *** 

c 

c SUBPROGRAM FCN REQUIRED eY THE ROUTINE ZSCNT 

C 

C at at at at at aj 4** afc^afra&atattatcafc^ajtattattafrafcaScafr******#**^ at*at * at * at * * *at at at at ** at at* 

SUBROUTINE FCN(X,F,N,T1 
INTEGER N , I , J ,NN 
REAL * 8 > ,F,R,U,A,B ,C,D,T 
CIMENSION X ( 4 1 , F ( 4 ) • T (41 
COMMON /FEFLEN/ R(41 

CCMMON /PARAM/ A ( 2 ,2 1 , B < 2 ,2 1 , C ( 2 , 2 1 , C (2 , 2 1 , U ( A 1 
Fill =-X ( 11*(U(11 +A (1 , 11 *X (3 ) *A( 1 ,21*X(A1 ) W R ( 11 
£-6(1 ,1)*M31-E(1 ,2l*X(Al )*T( 11 
F(2J=-X(21*(U(21+A(2,1)*X(31+A(2»2)*X(4))+(R(21 
£-B(2,ll*M31-E(2,21*X4All*T(21 
F(3) S -X(2)*(U(3J+C(1»11*X(1)+C(2»1)*X(2)1+(R(31 
£-D (1 ,11*X(11-0(2,11*X(211*T(31 
F ( Al -—X (Al*(U(Al+C(l ,21*X(11+C(2,21*X(21 1 +( R ( A) 

£-0(1 ,21* M 11-0(2 ,21* X (21 1*T( 41 
FETURN 
END 
C 

£***** ***** $$ at * * $ * *** * ** ** ************ * ** ** *** * *** * at ** *** 

c 

C SUBROUTINE TO TREAT SINGULAR CASE . WHEN EITHER 

C THE ••REVERSE" FLAG. REFUG CR THE ERROR FLAG, 

C ERFLAG IS SET , THE PARTiTiON INTERVAL, TS RILL BE 

C INCREASED. IN ADDI T ION .IF REFLAG IS SET, THE FAST 

C CHANGING COMPONENT WHICH HAS BEEN FOUND NOT TO 

C FLIP CORRECTLY WILL HAVE ITS SIGN CHANGED. 

C 

C* ** ** ***at * *** at a* * * * ** * **** * ** ** **************** a{ * Jfc * * * *** 

c 

SUBROUTINE TREAT! REFLAG, ERFL AG , CHf LAG 1 
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c 

c 

c 

10 



c 



1011 



c 

c 

c 

11 



c 

1012 

c 

c 

c 

12 



C 

1013 

20 



REAL * € XT, X , T,T S ,GMAX,XFTEMP ,SUM 
I MEG E F K , N, COUNT , CON , FAS T , QC CUNT , RGN , SI GN , I ER 
COMMON /TFEPAR/ XT (4 J ,X ( 4i , T ( 4) , TS,F AST , C ON , COUNT ,N 
6 »QMA X i QC OUNT , RON ,XFT EMP , SUM, SIGN , IER 
DO 10 S K = 1 ,N 

If ( .NOT.I ( (REFLAG.NE. 1) .AND . ( E RFLA G .NE. 1)1 
G.QR. (CHFLAG.EQ.l U J GC TG 10 
T(KJ=T (Ki+TS 
XI ( Ki - X ( Ki 

IF NEED TC SET X (FAST J=-XT( F AST ) S ERFLAG IS SET 

IF(.NOT.((REFLAG.EC.l) . AND . ( ERF LAG . EQ. 1 i JJ 
L GO TC II 

> (f AS7)=-XT( FASTI 
XTIFAST) =X(F AST) 

If {•NOT. (CHFLAG.EQ.O)J GO TO 1011 
T (K ) =T (K )+TS 
T (2i =T ( I J 
T (3i =T(2 i 
T(4i =T (3 ) 

T S=2 . 000 * TS 
CHFL AG=1 
T (K ) =T (KJ+TS 

CON= IDINT ( ( l.ODO-TI 1 i )/T S)+3 
CflUN T=1 
CCNTINUE 
P Ef LA G=0 
E RFL AG=0 
GO TO 2C 



JF ONLY ERFLAG I S SE T 

IF ( .NOT.! ERFLAC .EQ.UJ GO TO 12 
If ( .NOT. (ChF LAG.EQ.OJ J GOTO 1012 
Ti (K) =T (K ) +T S 
T 42 i =T ( 1 ) 

T (3 ) =T ( 2 ) 

T ( 4 J =T (3 ) 

TS= 2 . 0 DO*TS 
CFFL AG =1 
T(KJ =T (K J+TS 

CON= IDINT U l.ODO-TI 1 JJ/TS J +3 
COUN T = 1 
CCNTINUE 
PEf LAC =0 
E F.FL AG =0 



IF ONLY REFLAG I S SE T 
GO TO 20 

IF( .NQT.IREFLAG.EQ.1J i GO TO 20 
X (FA ST i=-XT (FAST i 

IFt.NOT. ( CHFLAG. EQ.O J 1 GO TO 1013 
T ( KJ =Tl KI+TS 
T ( 2J = T ( 1 J 
T(3J = 1 ( 2 ) 

T ( 4) = T ( 3 J 
T S=2 • ODOMS 
CHFL AG=1 
T (KJ =T ( K ) +TS 

CON= IDI N'T ( ( 1 .ODO-T ( 1 ) J/TS J+3 
COUN T = 1 
CONT 1NUE 
REFL AG=0 
EPFL AG=0 
CONTINUE 
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109 CONTINUE 

RETURN 
END 
C 

C*#** **$;£$ $*##*:$ *#*.*$ $.$$*:£ ♦*♦***# ****#$*£=*$******$**;<»$*$* 

c 

C SUBROUTINE TO DETECT FAST CHANGING COMPONENT OF 

C X.UPON RETURNINGTC MAIN PROGRAM, THE FLAGS ERFLAG 

C AND REFLAGWILL BE SET IF ANY SUCH COMPONENTS 

C DETECTEC 

C 

C 



SUBROLTINE DE TECT ( CX , R EFL AG , E RFLAG I 
REAL * 8 OMAX , X FT EMP , SUM , XT , X , CX ( 4),TS,T 
INTEGER K , CCG UNT, F AST , RON , SI G N , REF LAG , IER, ERFLAG 
£ , N, COUNT ,CCN 

COMMON /TREPAR/ X T (4 ) ,X< 4 J ,T ( 4 J , T S , FAST , CON, COUNT 
£ , N,OMA>,CCUUNT ,RON,XFT EMP , SUM, SIGN, IER 

C 

C IF DX (K J > OMAX FOR 10 TS ,THATCOMPON ENT IS CHANGING 

C FAST • 

C 

DO 105 K= 1 »N 

1F( .NOT. (DX< KJ.GE.OMAXJ) GO TO 1014 
IF(.NOT.( QC0UNT.GE.10 J J GO TO 1015 
F A ST=K 

XF T EMP=X( FAST I 
R0N = 2 
QC CUNT = 0 

1015 CONTINUE 

CC CUN 7=QC OUNT + 1 
1014 CONTINUE 

105 CONTINUE 

C 

C IF X ( FAST J=-X7( FAST) i*hEN FLIPPING, SET REFLAG 

C 

I F(.NCT. (RCN. NE.l ) ) GC TO 101o 
SUM-X1'(FAST ) +X ( FA ST ) 

IF I .NOT •( ( CABS I X ( FAST ) ) .LE . (0.2 CO* DABS 
(XT(FAST)J ) ) .OR. ( IER.GT. 0 ) )) GC TO 1017 
IF ( .NOT .( (DABS! X (FAST)) .LE. (C.3DC*DA£S 
( XT ( FAST) )) ) .AND. ( IER.GT .0) ) ) GO TO =0 
SIGN = IDI NT (XT (FA ST)/ X(FAST) ) 

IF( . NOT. (SIGN. LT .0) ) GO TG 1021 
REFL AG= 1 
CONTINUE 
I ER= 0 
ERF LAG-1 
GC TO 4 0 

IF (.NOT. (DABS(XiFAST)) . LE. ( 0 .3D0*DABS 
( XT ( FAS T ) )) ) )GG TO 31 

S IGN=IDINT( XT (FAST) /X (FAST) ) 

1F(. NOT. (SIGN.LT .0) ) GO TO 10 2 2 
R EFLA G-l 
cqnt INUE 
GO TO 40 

1 F ( .NOT .( IER.GT. 0) ) GO TO 40 
I ER= 0 
ERFL AG=1 
CONTINU E 
CONTINUE 
10 lu CONTINUE 

C 

C IF NO COMPONENT IS CHANGING FAST, BUT THERE IS ERROR 

C FRCM 2SCNT 

C 



£ 

£ 

1021 

30 

£ 

10 22 

31 

40 

1017 
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I F ( • NOT « ( (RON . EQ. 1 J . AND, ( IER.GT.OJJJ GO TO 1020 
I ER = 0 
EFFLAG=1 

10 20 CONTINUE 

RETURN 
ENO 
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APPENDIX D 

DERIVATIONS OF THE RELATIONS BETWEEN £ , £, AND STABILITY 

* y 



to 



Neutral Stability 

In this case, z - z =0 and equation 4.5 reduces 

x y a 



n l (x ' x el } + + y l Cy '>'ei :i = 0 



It follows that both sets of hyperbolas merge into one and 
all the points on the common hyperbola are equilibrium 
points. The lefthand side of the second condition in equa- 
tion 4.4 can be manipulated as follows 



(u + ay 
= ac 
= ac 



e )(v+cx e ) - (b+ax 



(n 1 + y e )(M 2 + x e ) 



x e^ n r n 2^ + 



e VM 2 



e )(d-cy e ) 

(vV x e )(n 2 + >’e } . 

- yi ) + n l( J 2 - r ' 2 w 1 J 



= 0 



This result implies that the constant term of the character- 
istics polynomial D(s), is zero ; hence one of the eigenva- 
lues, s equals to zero. Factoring out the characteristic 
polynomial, we get the other eigenvalue as 



S 



2 



( U +ay e ) + (v+cx e ) 
a (ni + y e ) + c(y 2 +x e ) 



(eqn D.l) 
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For points on the first quadrant hyperbola, x > - n and y 
> -r\ ; therefore s 2 <0 and neutral stability exists. 

Conversely, s 0 > 0 on the third quadrant hyperbola which is 

therefore unstable. The results are summarised as fellows : 

(1) Khen £ x ~ z y ~ 0/ inf initely-many equilibria 
exist as points cn the two hyperbolas on which x = y = 
C; 

(2) The first guadrant hyperbola is neutrally stable; 

(3) The third guadrant hyperbola is unstable. 

2. Inters ec tion s in First and Third Qu adra nt 

The proofs for the following results are given in 
this section : 

(1) when both equilibrium points are on the first 
guadrant hyperbolas, one is stable and the other 
unstable 

(2) when one equilibrium point is on the first guad- 
rant hyperbola and the other on the third, both can be 
unstable or one w ill be stable and the o th er unstable. 

The straight lines given by equation 4.7 are plotted 
on the e x , Zy plane as shown in Figure D.l. It also shews 
the corresponding regions on the e x , e y plane as x e2 and y 
vary. 

let (x , y el ) be the first equilibrium point on the 
first guadrant hyperbola. The first stability criterion in 
equation 4.4 is automatically satisfied since 

( u+ ay e ) + (v+cx e ) 

= a ^i + X e ) + c(u 2 +x e ) 

> 0 for * e2 > -u 2 and y e2 > -n : 
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Figure D.1 Effect of Varying x „ and y on e , e Plane. 

e 2 e 2 x y 

Conse guently , the stability of (x el , y^') is solely deter- 
mined by the second condition in equation 4.4. The second 
condition can be rewritten in the following manner : 



(u+ay )(v+cx ) - (b+ax )(d+cy ) 
J e e e e 



= ac 
= ac 
> 0 



*e (n i” n 2 ) + y e (y 2” 1 J l^ + ' J 1 (n 1 -n 2 ) + n i (y 2“ y l^] 

■ E y lx e t,, l ) + 



or 



s 

£ y 5 
7 u 



£ 

X 



<yr e > 

( y x e ) 



(eqn D . 2) 
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9 



For (x 



el 



line of stabilit 
+ u x ) and (y el + 



Y el ) r e 3 ua tion D. 2 represents the 
and this line has a slope between 
ri l)/y 1 . 1 1 is shown in Figure D.2. 



boundary 



n x /(x 



el 




Figure D.2 Boundary line of Stability on e Plane. 

To investigate the stability of (x^, y e2 ), we 

substitute them into equation D.2 and obtain the condition 
for stability as 



u 

e y J e x 
s 



(Pl + x el ) 



(eqn D.3) 



100 



Equations D.2 and D.3 are only different in the 
directions cf inequality. This means that on one side of the 
line e y = £ x ( n 1 + Y e l) / ( ,J 1 + x el ) , one e i u ilibr iu m point is 
stable and the other is unstable. On the other side of the 
line, the opposite conditions exists. 

Note that (x e2 , Y e 2^ ma Y or ma Y not satisfy the 
first condition of ecuation 4.4. This condition is 

a (ni+y e2 ) + c ^! + e x +x e2^ > 0 ^ eqn D - 4 ^ 



From Figure D.1, equation D. 4 and earlier results the 
following deductions can be made 



(1) In the first quadrant of the e x , 



plane, e > 0, 



*.2 > "V 



y e2 > ; 



hence equation D.4 is satisfied. 



The region labelled 

(x 



m. 






in Figure D.3 is stable for 
e2' 1 e2 ) iUt Unstabl e fcr ( x e i' Y e 1 ) ; 

(2) In the third quadrant of the e x , plane but 

between the two lines where x . > 0 and y _ > -n. , 

equation D.4 is again satisfied. The region labelled 
k||| is also stable for (X , Y e2 ) but unstable for 
) ; 

e plane, 

- y r 

e 






(x 



el' J el 

<3) In the second quadrant of the £ 

v _ , 1 < 0, 

1 e2 1 x 

fied ; so both (x , , y .) and (x 0 , y «) are unstable. 

6 1 6 1 0 Z 0 Z 

The region is labelled 



x ~ < 

x y * e2 

equation D. 4 is not satis- 




in Figure D.3; 

(4) In the fourth quadrant of the e x , plane, equa- 
tion D.3 is not satisfied ; so (x „ , y „) is unstable 

0 Z J 0Z 

but (x el , y el ) is stable. This region is labelled 
in Figure D.3; 

(5) In the region labelled 




e quat ion D.3 is 



not satisfied ; so (x e2 , y e2 ) I s unsta kle tut (x e j. • 

y ) is stable. 

J el' 
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^ x el ’ y el^ 
unstable; 
Cx e2> y e2 ) 



stable 
as above 



both 

unstable 




(x el’ y e p 
stable; 

(x e2> 

unstable 




as above 




x 





(y l +x el 



) 

) 



Figure D.3 Regions in e , e Plane- 

x y 

By carefully noting the signs of (x e2 , y e2 ) and the 
various regions in Figure D.3, all the previous deductions 
can be combined ; we conclude that : 

(D When both equilibrium points are on the first 
guadrant hyperbola, one is stable and the ether 
unstable; 

(2) When one equilibrium point is on the first guad- 
rant hyperbola and the other on the third, both can be 
unstable or one will be stable and the other unstable. 
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3 • T wo Equilibria in the Third Quadra nt 

The simplest way to show that both equilibrium 
points in the third quadrant are unstable is t o refer to 
Figure D.4. Clearly, we must have x < - v/c and y < - u/a 
for both equilibrium points. In that case, the first 

stability criterion in equation 4.4 is not satisfied since 

(u+ay ) + (v+cx ) < 0 

0 O 

Therefore, both equilibria are unstable. 



-x _ v. b. 




Figure D. 4 Two equilibria in the Third Quadrant. 



4 . Rep eat ed Equil ibri a 

This case corresponds to operating exactly on the e 
= e (y + n,)/U + y J on the e , e plane. The proof is 

obtained by substituting x^ = x n , y „ = y , into equation 

u 3 e2 el i e2 ‘el 

4.6 and solving for £ y . 
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(eqn D.5) 



y 

y e 2 = y el ■ ^ 


( x e i*Pi) 


' n i 


£ 


X 

X 0 = X . = — 

e2 el e 


( y el +n l ) 





Rewriting eguaticn D.5, we have 

£ x y el = £ y x el + 1 , l c y ' n l e x 
£ y x el = e x y el + n l £ x ' v l € y 



subtracting one from the other. 



2e y = 2e x , + 2y,e - 2n-,e. 
x'el y el 1 y i - 



e 



y 



(y e i + V 



( X el + V 



(eqn D.6) 



Thus we have shown that the case of repeated equi- 
libria corresponds to points on the straight line indicated 
in Figure D.3. 

Equation D.6 can be also obtained by setting the 
second stability criterion egual to zero ( see equation D.2 
). That is equilavent to saying one of the eigenvalues, s^ 
is equal to zero. The sign of s 2 is determined by consid- 
ering the first stability criterion. Following the same 
argument as in neutrally stable case, we can prove that 
repeated equilibria on the first quadrant hyperbolas are 
neutrally stable. On the other hand, repeated equilibria on 
the third guadrant hyperbolas are unstable. 
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APPENDIX E 



RELATION BETWEEN STABILITY AND e , e PLANE: VERIFICATIONS 

x y 

Seme representative points on the e x * plane are 

chosen and the corresponding stablities of the equilibria 
calculated. The points chosen are marked F, G, H, M, N , P, 
Q, T and W in Figure E.l. 




Figure E. 1 Experimental Verifications. 
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the £ 



we 



Fr oced ure 



Having chosen the points on 
have to work backward to obtain a. 
Suitable (x el , y el ) are then chosen, 
tion cf r, s, X g2 and y g2 . (x g2 , 
directly from eguation 4.6. From all 
eigenvalues can be calculated and 
theory. 



X , £ y plane, 

b, c, d, u and v. 

followed by a calcula- 

y i can be obtained 
e 2 

these parameters, the 
results compared with 



2. Fesults 



The results are tabulated in 
with the theoretical results given in 



Table VII. 
Figure E. 1 . 



They agree 
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TABLE VII 

Results of Experimental Verification 
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4.2(d) 



APPENDIX F 

PROGRAM TO PLOT TRAJECTORIES IN 1*1 SYSTEMS 



C 
C 
C 
C 
C 
C 
C 

£********* *** ******** ******** ********* *********** 

C PROGRAM 1 1PL FORTRAN 

C THIS PROGRAM RLOTS THE TRAJECTORIES OF A 1*1 SYSTEM. 

C IT CALLS AN IMSL ROU TIN E , DVERK FOR INTEGRATION. 

C EEFOFE EXECUTION , DETERMINE THE SCALE OF THE 

C PLOT AND NUMBER OF CURVES REGUI RED. 

C THE VALUE , "STEP" IS DETERMINED EY THE RANGE OVER 

C WH IC H HE V 1 A NT TO SET THE INITIAL POINTS. CHECK A,B, 

C C»D» ANO U ( 2 ) • FI X EITHER X(1J OR X ( 2 i AND VARY THE 

C OTHERJTHE ONE FIXED HAS TO BE RESET AFTER EACH 

C CURVE HAS BEEN DRAWN. IT IS RESET IN THE LAST 

C ASSIGNMENT STATEMENTIN THE DO LOOP. TO EXECUTE, 

C ENTER "IlFL". THE EXEC FILE,"I1PL EXEC" 

C MUST EE IN THE DISK. IT MUST BE EXECUTED ON A 

C TERMINAL ATTACHED TO THE TEKTRONIX 616. 

Q***** ******** * * ******************************* *********** 

c 

c 

c***** **************************************** ************ 

C VARIABLE CEFINITIONS 

£***** * *** * **** 4 ***^** ************************ * *********** 

C >= VE CTOR OF LENGTH 2 CONTAINING ThE UNSTABLE 

C EQUILIBRIUM POINT. 

C AfB,C,D,U,R = ATTRIT10N CCEFFI ENT S • 

C CAPX tCAPY fCAPlX, CAP1Y- ARRAYS USED TO STORE X"S 

C FOR PLOTTING PURPOSES. 

C CEL= INTEGRATION STEP SIZE. 

C IND, CC, W,TOL,TEND=PARAMETRS REQUIRED BY IMSL ROUTINE 

C OVERK 

C IER= ERROR MESSAGE NUMBER FRUM DVERK 

C 01 R= CONSTANT FACTOR DETERMINING THE DIRECTION 

C OF PERTURBATIONS FROM THE UNSTABLE POINT 

C FCN1 = EX T ERNAL FUNCTION REQUIRED EY DVERK 

C NSTE P = NUM BER CF CURVES TC PLCTjALSO s OF STARTING 

C STEP=INT EFVAL BETWEEN DIFFERENT1 N1TI AL POINTS. 

C INTEGRATION POINTS. 

C K , KM 1=C 0UN7 ER S FOR PLOTTING ROUTINE 

Q***** ******** ************ ************ ********* *********** 

c 

INTEGER N,IND,NW , I ER , K , NPOI NT ,KM 1 , KK , J J , K 1 , K2 ,KM2 

INTEGER N , INC ,NW , I ER ,K, NPOI NT , KM1 ,KK , J , NSTE P 

PEAL * 4 X(2),CC(24J ,W(2»9),T»T0L, TEND, DEL, R(2),A,B» 

£ C , C» U 12 J ,CAPX(1000J, CAP Y( 10001 , X A < 10 00 i , X El 1 0 CO I , STE P 
EXTERNAL FCN1 
COMMON R,A,e,C,D,U 

OATA NW/2/ ,N/2/, T/O. 0/,T0L/0 .001 0/ , I ND/1/ ,IER/0/ 

C 

C ENTER ATTRITION COEFFICIENTS HERE 

C 

A= C. 6000 

e=o. sooo 

C=1.0000 
0=5. CC 
U( li =0.60 
U(2) = 1.00000 

c 

C ENTER ONE EQUILIBRIUM POINT HERE FOR THE CALCLLAT ION 

C OF R VECTOR 

C 

X( 1J=5.0CCOO 
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X 

p 

p 



2) =4. C C C CC 

1J=X(1 J*(U(1J«-A*X( 2))*6*X(2J 
21 =M2)*IU12J + C*X( 1 )J*C*X( II 



C 

C ENTER R DIRECTLY HERE IF THEY ARE KNONN;R£MGVE 

C THE COMMENT CHARACTERS, 

C 

C p(l]=444** 

C P(2)-***** 

C 

4 444 44*44 4 444 44* 4444 444444444444 44444444 44444 4444444 

C SET PLOTTING PARAMETERS 

C44444 4 4 44 4 444 4 4 4 44 44 4 4444 4 44 4 44444444 4444 4 444 4 444 4 4 44 4444 

c 

CALL TE K £ 1 £ 

C CALL COM P RS 

C CALL VRSTEC (0,0, 0) 

CALL PAGE (14.7, 11. 5J 
CALL NQEPCR 
CALL CROSS 
CALL BLOWUP (0.5) 

CALL AREA 20 (12. 7 ,9.50) 

CALL XN A ME ( ' TCTA L X FGRCES',100) 

CALL YNANE(' TOTAL Y FORCES', 100) 

CALL FRAME 
C 

C INSERT HEADINGS WITHIN CLOT ES ;R EMOVE COMMENT 

C CHARACTER 

C 

C CALL HEAO IN (' 444444444444444444444444 4444 444$ • 

C £,100,1. ,4) 

C CALL HE AD IN ( ' 4444444 4444444444444444444444$ • 

C £,100,1. ,4) 

Q CALL HE A D IN ( • *** 44 44 44444 4444 44 4 44$' , 100 

C £.100,1. ,4) 

C CALL HE A 0 I N (' 44 4 4 4 4 4 4 4 4 4 4 4 4 4 444 4 4 4 4 4 4 4 4 4 4 4$ 1 

C £ ,100,1., 4 j 

CALL GR A F (-3. 000 , 4 SC ALE' ,10.0,-6.00, 4 SCAL E' , 1C.00) 



C INSERT INITIAL CONDITIONS FOR INTEGRATION 

C 

C 



C 

c 

999 



STEP = 0. 5CCC 
NSTE P = 1 0 
X( 1) *1.0CCO 
DO 100 J- 1 ,NS7EP 

X (21 =-5.00C+FLQAT ( JJ*STEP 
K - 1 

DEL=C.01C 
NPOI NT = 200 

I F(.NQT. ( < IER.LE. 0 J .AND. ( IND . GE.O ) .ANC.t K. LE. 
i NPOINT))) GO TO 6 

TEND=FLCAT (K)*DEL , 

CALL OVERK (N,FCN1 , T ,X, TEND , TOL , IND , CC , N W ,W , IER i 
X A ( K J = X ( 1 J 
XB ( K ) = X (2 ) 

WRITE! 6 ,99 9 J IE R 



,3X, 
1 J 
21 



FOPNAT ( 'O' 

C A P X ( K ) =X ( 
CAFY(K)=X( 

K = M1 
GO TO 5 
CONTINUE 
K M1 = K- 1 

CALL CLRVEICAPX.CAPY 



IER= 4 , I 3 J 



KM1 . 3 i 



C 



on 



RESET INTEGRATION PARAMETERS AFTER EACH CLRVE. 



I E P= 0 
I ND= 1 

X < 1J = 1.0C0C 
100 CONTINUE 

CALL END F L ( 0 J 
CALL CONE PL 
STOP 
END 
C 

C SUBROUTINE FCN1 REQUIRED BY OVERK 

C*** 3 ^* ** ** *$$$< <$$#<*****# $*£#**** * *** $ **♦■<:** 

c 

SUBROUTINE FCNH N ,T, X ,XPRIJ 
INTEGER N»lfJ 

REAL * 4 X(N) .XPRI (NJ ,T,R (2) , A ,B » C ,D, U (2 J 
COMMON ft ,A ,B, C,C, U 

XPFI (1 )=-X(l J *(U( 1) + A*X(2JJ^R(1J-E*X( 2) 

XFRI ( 2 i =-X (2 i *(U( 2)+C*X(l JJ+R(2J-C*X( 1J 
RETURN 

END 
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APPENDIX G 

PROGRAM TO PLOT BOUNDARY CURVE 



C 
C 

c 
c 
c 
c 

C***** $****$*$ i *4 $$ $$*#*$$ $$*:{: $*4:$ $**4: ***$*$ *$****$.$.**# 

C PROGRAM EK1 FORTRAN 

C THIS FRO G F AM DETERMINES THE EOUNDARY SEPERATING THE 

C COMA INS OF ATTRACTION IN A 1*1 P RCBL EM . B ACKW ARD 

C INTEGRATION (NEGATIVE TIMEJ IS USED STARTING 

C FROM AN UNSTABLE EQUILIBRIUM POINT, 

C BEFORE EXECUTION CHECK A , B, C ♦ D , U ( 2 1 , AND X(2). 

C TO EXECUTE ENTER BK1.THE EXEC FILE,"6K1 EXEC’ 1 

C MUST EE IN THE DISK. IT MUST EE EXECUTED ON A 

C TERMINAL ATTACHED TO THE TEKTRONIX 618. 

C***** **** **** * * * * £$$$ 4> ******** ***** **** **** ******* 

C 

c 

£*****4***4‘******** 4** * *** ***************** * ** * **** ******* 

C VARIABLE DEFINITIONS 

£***** ******* '********* **** *** ***** **** ********* *** * 4 ****** 

C X=VE CTOR CF LENGTH 2 CONTAINING ThE UNSTAELE 

C EQUI LIBRIUM POINT . 

C A«8»CjDtU»R = ATTRITION CCEEFI ENT S . 

C CAPX ,CAPY,CAP1X, CAP1 Y=ARRAYS USED TO STORE X"S 

C FOR PLOTTING FURPOSES. 

C CEL- INTEGRATION STEP SIZE. 

C JNDt CC» W # 7GL » TEN D=PARAMETRS REQUIRED BY IMSL ROUTINE 

C DVERK 

C IER= ERROR MESSAGE NUMBER FROM DVERK 

C DIR= CONST ANT FACTOR DETERMINING THE CIRECTION 

C OF PERTURBATIONS FROM THE UNSTABLE POINT 

C FCM =EX T ERN AL FUNCTION REQUIRED BY DVERK 

C E=CONST ANTS EQUAL l.OiFOR SETTING DIRECTION OF 

C CF P ERTUR BAT I CN TOGETHER WITH DIR 

C- K • K1 *K2= COUNT ERS FOR PLOTTING ROUTINE 

£********* **** * ****** ******** ********* ******* ************* 

c 

INTEGER N,IN0,NW . I ER ,K, NPOI N T , KM 1 , KK ,JJ,K1,K2 ,KM2 
REAL * 4 X(2J ,CC (24) ,W(2,9J , T ,TOL, TEND ,D EL, E ( 2J , R, A, B 
£*C»D»U»CAF>(1 COO ltCAPY( 1000 i , 

GCAPlXtlOOCl tCAPl YdOOCJ ,CIR 
EXTERNAL FCN1 
COMMON R ( 2 J 

COMMON / PARA/ A,B*C»C,U(2) 

OATA NW/2/ fN/2/#T/0.0/tT0L/0.0010/»IND/l/t I ER/O/ 

£ ,E/+1. 000 ,41.000/ 

A=0. 7 
E = 0. A 
C = l. 0 
0=0. 6 
U( 1J=0. 15 
U(2J =0.2 
C 

£***** * ******* 4 * ****************************** ************ 

C CALCULATE R,S VECTOR USING DATA FROM LP PROGRAM 

£************* 4 ******* ************************************ 

c 

X( 1) =£.70 
X ( 2J =2.00 

R ( 1) = X ( 1 i*(U(l) + A*X( 2JJ+E*X( 2) 

R(2J=X(2 J«(U(2)+C*X( 1JJ+0*X( 1) 

WRIT E ( 6 ,9 998 HI , X( U , I, R( II , 1=1, Ni 
9998 FORMA7( , 0*»3X ? , X(* » I 3 #• J=*»F12.4»4X» , P( I t 13 » 1 )=* ,F12.4 

C 

£***** ******** ************************************** 4 ***** 



111 



C INSERT INITIAL CONDITIONS FOR INTEGRATIUN 

C$ 4‘4< ** 4 4** 44'4‘444444 44‘4s44‘44444444*4‘4444444«4‘444<4‘4444‘4‘ 4 44 4*4*4 

c 

Kl = l 
K2 = I 

CO 100 J J= 1 » 2 
C 

C FOR J J - 1 £ 2 » BACK W ARC INTEGRATE ALONG ERANCH * 1,2 

C RESPECTIVELY 

C 

I F I .NOT . I J J. E Q. 1) ) GO TO 50 
0 1 P= 0 . 03 

X( 1 JsX(lJ + DIR*E(l ) 

XI 2J = X(2J+D1R*E(2J 
GO TO 60 

50 CONTINUE 

D I R=-0 . C3 

XI 1 1 =6. 70* D 1 R* E 1 1 J 
XI 2)=3.000+DlR*EI2) 

60 CONTINUE 

K = 1 
T = 0 . 0 

0 EL=Q *100 

IF (JJ.EC. 2) D EL=0 . 100 
NPOI NT=200 

1 ND= 1 
C 

C4c*444< 44444>4444 4 4 44 44 4 44 44 444444444444 $ 44^44444 444 4 444444 4 

C 

C WHILE NO ERROR AN C NUMBER OF INTEGRATION STEPS 

C LESS THAN NPOINT 

C 

C44444 *$*-***$4' 44444 44 444 44 4 44 4 4444444 4= 4* 44 * 4 444 * 444 4 4444444 

c 

5 I FI. NOT. I I IER .LE. Cl.ANU. I IND. GE.O J .AND. I K. LE. 

£ NPOINT ) J 1 GO TO 6 

TE ND=- I FLO AT IK J *0 EL 1 

CALL INTEGRATION ROUTINE 

CALL 0VERK(N,FCN1,T,X, TEND , TOL , IND , CC , N W, W , I ER 1 

PRINT ERROR MESSAGE IN ANY 

WP ITE<6,999) IEP 
FORMAT I *0 * , 3X i • IER** ,13} 

STORE XU) » X 1 2 ) FCR THE TWO BRANCHES OF THE 
BOUNDARY IN TWO ARRAYS 

IF I.N0T.IJJ.EQ.1J )G0 TO 501 
CAPXIK) =X(1 J 
CAPY IK I -X 12 ) 

K 1 = K1 + 1 
GO TO 601 
CONTINUE 

CAP1XIK ) = XI 1 J 
C AP1YIK J^XI 21 
K 2 = K 2+ 1 
CONTINUE 
K*M 1 
GO TO 5 
CONTINUE 
CONTINUE 

SET NUMBER OF PLOTTED POINTS TO ONE LESS THAN 
THE NUMBER OF DATA POINTS 



£ 



C 

C 

c 

999 

C 



501 



601 



6 

100 

C 

C 

C 
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OOOOOOOOOOO oo 



KM 1= K 1— 1 
KM2=K2-1 



C***c ****** **** * 

C COMMENCE 



*********** * ******* *** ********* * **** ******* 
PLOTTING THE BOUNDARY 
C***** * *** **** ***** ************************ *** * *********** 
C 



CALL TEK6IE 
CALL COM F F S 
CALL VR S T E C ( 0 ,0 , 
CALL FACE (14.7,1 
CALL NOB R C F 
CALL CROSS 
CALL BLOWUP (0 .55 i 
CALL AREA2D(13.7, 
CALL XN AM E ( ' 

CALL YNANE < • 



0) 

1-5) 



9.0) 



TOTAL X FORCES ' .100 ) 
TOTAL Y FORCES' ,100) 



INSERT HEADINGS WITHIN QUOT ES ;R E MGV E COMMENT CHARACTER 



CALL 
E ,100 
CALL 
& ,100 
CALL 

& y 100 

CALL 
E ,100 
CALL 
CALL 
CALL 
CALL 
CALL 
STOP 
END 



FEA 
, 1 • , 
HEA 
» 1 • » 
HEA 
» 1 « » 
HEA 
* 1- , 
GRA 
CUR 
CUR 
ENC 
DON 



DIN ( ' ******* ******** *** ******* ****** $ 1 
A) 

Q ] N( 1 ******* ************ **********$* 

A ) 

0 IN ( ' *** **** ************ **$' ,100 
4 ) 

D 1 N i * *********** ************* *****$* 

4 ) 

f <6.000, 'SCALE' ,7.000,0. 400, 'SCALE* ,4.000) 
VEiCAPX, CAPY , KM 1 , 2) 

V E (CAF1X ,CAP 1Y,KM2, 2) 

FL (0) 

EFL 



£***** ** ****** 
C SUBPROGR 

C***** ******** 

c 

SUBRQUTI 
INTEG 
R EAL 
COMMO 
C CMMO 
X PRI ( 
X PRI i 
RETUR 

END 



* 


* 


* ***** 


** ** 


* 


** 


* 


*** 


A 


M 


CALL 


BY I 


N 


TE 


G 


RAT 


* 


* 


****** 


**** 


* 


** 


* 


*** 


N 


E 


FCN1 l 


N , T, 


X 


,x 


PRI) 


E 


R 


N , 1 , J 












* 




4 X (N } 


, X PR 


I 


(N 


) 


,T, 


N 




R ( 2 ) 












N 




/PARA/ 


A ,B 


, 


c, 


D 


, U( 


1 


) 


= -X(l) 


*(U( 


1 


) + 


A 


*Xl 


2 


) 


= -X(2) 


*( U( 


2 


)* 


C*X< 



N 
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c 

c 

c 

c 

c 

c 

c*** 44 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

£44444 

c 

c 

£44444 

c 

£♦ 4444 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



PROGRAM TO 
OPTIMUM 



APP E NO 1 X H 

0 QTA11A PAYOFF MATRIX 
MIXED STRATEGIES 



AND 



4444 444444*44 $4 4444444444 444444 jjt 

THIS PROGRAM COMPUTES THE PAY 
OPTIMUM MIXED STRATEGIES FOR 
IS BASED ON THE THEORY GIVEN 
THE P A 't OFF FUNCTIONS ARE COMP 
BUT ONLY 15*15 SAMPLED MATRIX 
THE OUTPUT ALSO I NCLUCES MATR 
F 1NT 1 M . 

BY CHANGING ALAMDA TO LESS TH 
OBJECTIVE FUNCTION CAN BE MAD 
MORE OF PINT I M IN ACCORDANCE 
A ( X » Y I =ALA MDA* ILY-LX)+(1-A 
BEFORE EXECUT ION,ChECK ATTRIT 
QX,QY ,AL AM CA, X1S,X2S, STEP 1, ST 
LABELLED C1,C2 AND SCALE IN T 
ROUTINE. PERTURBATIONS ARE GIV 
LABELLED C3 AND CA. 

TO EXECUTE, YOU NEED THE EXEC 
AND A TERMINAL CONNECTED hi TH 
11 LGSQP" WHEN READY. 

4 4 44 4 * * 4 4 4 4 4 4 444 4444 444444444444 



44444444444444444444 
CFF MATRIX AND 
X AND Y. THE ALGORITHM 
IN CHAPTER 5. 

UTED FCR 60*60 X,Y 
IS PRINTED. 

IX FOR L X ,LY , AND 

AN 1 , T H E 

E TO EMPHASIZE 

WITH: 

LAMD A J *F I NT I M 
ION COEFFIENTS, 

EP2, CARDS 
HE PLOTTING 
EN BY CARDS 

FILE "LOSCP EXEC 11 
A TEK616. ENTER 

4444444444444 4444444 



4444 4 444 4 4 4 44 444 44 44 4 44 4 4 44 4 4444 44****** 

VARIAELE DEFINITICNS 

4444444444444444444444444444444444444444 

X ( 2J = A PR AY CONTAINING X AND Y 
XIS, X2S=VARIBLES USED TO SET THE FART 
OF X AND Y VALUES TO CALCULAT 
A(X,Y). THERE ARE 60 INTERVALS 
Y IN THE RANGE 0.2Q AND 0.75Q 
S TEP 1 * S 7 EP2-T 0 GBTAIN THE CORRECT STE 
AS DESCRIBED ABOVE 

XFRIN1 ,XPR IN2 ^VARIABLES USED TO SAFER 

VALUES OF X AND Y 
TX,TY, R,S = A$ DEFINED IN CHAPTER 5 
RX,RY=(OX-X) , (QY-Yj RESPECTIVELY 
XL OSS »YLOSS=LX»LY IN CHAP. 5 
C C »NW »W,TQL»T END, DEL, IND* FARAMETE P S R 

IMSL ROUTINE, DV 
IED=ERROR CODE FR CM DVERK. 

ZMAT= MATRIX USED FCR PLOTTING SURFACE 
R ATL 0 S = F A YCFF MATRIX 
F INT I M= F IN ISH TIME MATRIX 
R ATA V E= 1 5* 15 PAYOFF MATRIX FOR OUTPUT 
X LGSP R = L X MATRIX FCR OUTPUT 
YLOSPF=LY MATRIX FCR CUfPUT 
F T PR= F INTI M MATRIX FCR OUTPUT 
ROWM 1 N = RGW MINIMUM IN PA Y CFF MATRIX 
COLMA X = COLUMN MAXIMUM IN PAYOFF MATRI 
AXMI N=MAXINUM OF THE THE ROW MINIMUM 
X INMA X = MIN IMU M OF THE COLUMN MAXIMUM 
RHS= ARRAY CONTAINING RHS CF CONSTRAIN 
0 ECO EF= A RRAY CONTAINING COEFFICIENTS 
FUNCTION IN LP. 

A LAMDA=NUMeER BETWEEN 0 AND 1 TO CETE 
RELATIVE EMPHASIS GF LOSS AND 
IN AIX.YJ 

S MMI N = S MA LLE S T VALUE OF RCW MINIMUM 
PAY=15U5 MATRIX USEC TO SAFEKEEP RAT 
CCMPUTE THE OPTIMIZED STAREGIES 



4 44444444444 
4 4444 4444444 

I TONS 
E 

FOR X AND 
P INTERVAL 
E EP INITIAL 



ECU I R ED BY 
ERK. 



A ( X , Y } • 



TS 

OF OBJECTIVE 

RHINE 

FINTIM 



AVE AND 
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nonooooooooooooo 



NS1,NS2=DIMENSI0N OF PAYOFF M ATR I X , EQ UAlS 60. 

N CGU = F LAG TO INDICATE END OF 1NTEGRAT I CN. 

JJ,KK= INDICES FOR MATRICES 
C QUNT = COUNTER FCR INTEGRATION 
KKM1=MJMBEF OF POINTS PLOTTED 

FCN1 ,2 ,3 ,4=SU BRGUTINES REQUIRED BY DV ERK ; TEE Y 

CONTAIN 1*1 EQUATIONS FCR TFE VARIOUS 
STAGES FOR THE BATTLE 

T =TI M E 

X A,YA ,ZA = ARRAYS USED FOR STORING AND PLOTTING 
XPRI N 1 » XPRIN2 »RATLOS AS DEFINED ABOVE 
Z COUN 1-CQUNT E R FOR X A ♦ YA 
J C = CO UNT ER FOR ZA 
T IKE 1 = T 1 IN CHAP* 5 
T 1 = T 2 IN CHAP. 5 

N PLOT = NUMB ER OF PLOTS REULlRED 
£***** **************** **** ************************* *** **** 
REAL X 15 ,X2S, X(2) ,R,S,RX,RY,TX,TY,CX,CY,A,C,U,V, 

£ W (2, 9 I ,T,TOL, TEND , CEL , XA ( 3650 ) , YA ( 365 C J , RAT 1 0 
£ , RATLCS (60 ,60 J .XPRIN1 ,XPK IN2 , XT1 , XT 2, RGTRAT 

£ , ZMAT UO ,60) , RORMINI 60 J ,C OLMA X ( 6 0 J ,AXN IN, X INMAX 

REAL F INAL .FI NT IM ( 60 ,60) , TIME 1 , RAT AVE ( 20,201 
£ , SMMI N,PAY (17 ,3 2) ,TEMPG,RhS(15) ,OBCOcF (15 ) , AL AMOA 

£ , XL0S(6C,6C) , Y LOS (60,60) , F TPR ( 20 , 20 i , 

£ Y LOS PR. (20 ,20 J ,CC( 24 J ,T1,STEP1 , ST EP2 , Y L CSS , XLOSS 

REAL YLCXL ,ZA ( 36501 ,XLGSPR(20 ,20 J 
INTEGER 1ND,NW,IEC,K, COUNT »NSTEP»NS1#NS2,J ,NOGO 
£ , ZCUUNT ,NX ,NY ,NPUINT , J J ,KK , JK , JC , KKM1 

COMMON /PAPA M/ A, C ,U , V ,S , R ,B , L 
EXTERNAL FCN1 
EXTERNAL FCN2 
EXTERNAL FCN3 
EXTERNAL FCN4 
C 

C*** *********** ************************ ******************** 
C SET PARAMETERS ANO COUNTERS 

C***** **** **** * < * * * *** **************** ******** * *********** 

c 

A= 0. 7 0 
C = 1.00 
U=0. 15 
V = C. 200 
B=0. 4000 
D = C. 60 0 
A LAMO A= 1. 00 
F INAL-O. 10 
N W=2 
N = 2 
T = 0. 0 
TCL= 0.001 
I NC= 1 
D E 1=0 .Cl 
C Q UN T = 1 
QX = 10 .CC 
Q Y = 7 • C C 
X 1 S=3 • C 
X 2 S= 3 . C 0 
JC=0 
NQGO=C 
N S 1=6 C 
N 52= 6 C 
N 5=N S 1 - 1 
ZCOUNT=Q 
NPOI NT=NS1*NS2 
C 

c***** ******** ***** ******* ************************** ****** 
C START COMPUTING 60*60 PAYOFF FUNCTIONS 
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£***** ******** * ******************************************* 

c 

DC 100 J - 1 »N S 1 

ST EFIM FLOAT (J-NS1 11/39.33 
DO 101 K= 1 t NS2 
JC=JC+1 
1 * 0.0 
COUN 1= 1 
IN D- 1 

STEP 2= ( F L CA 7 ( K - N S2 ) 1/54.428 

C CARD Cl: 

XI 11«X1S*(2 .333 + STEPli 

C CARD C 2 : 

Xl21=X2S*(l.S17+STEP21 
XFRIN1 = X( 11 

> PR IN2= X ( 21 

X A ( J Cl = XPRI N1 
Y A ( J Cl = XPRI N 2 
R = X( 11* ( U+A*X( 21 1+B*X( 2 1 
S=X ( 2) * { V+C*X( 111+D*X( 1 J 

rx^ox-x( u 

RY = GY-*X( 21 

C CARD C3 : 

> U1=XPRIN1-0.01 

C CARO C4 : 

X(21=XPRIN2-0.05 
T X=RX/R 
1 Y =R Y/S 
C 

Q**** ********* * ******* ************************ * *********** 

C STAGE 1 

C***=** ******** ******** ************************ ***** ******* 

C 

7END=0. 0 

1 1 ME 1=A Ml N1 ( TX » T Y J 

5000 1 F ( .NOT •( (TEND.LE.T IME1 i .AND. (NOGO. EC. 01 11 

£ GO TO 6000 

7END=FL0A7(CCUNT1*DEL 

CALL DVERK(N t FCN4,T, X ,TE ND ,TOL, I ND,CC 
6 t NW » W » I E D 1 

I F( .NOT. UX(1).LE.( FINAL* CXI 1 .OR • (X ( 21 
£ . LE . ( F IN AL*CY 1111 GO TO 1023 

I F ( .NOT. ( X < 1 1 .LE. (FlNAL*gXll 1 GC TO 52 
XLOSSM l.-FINALl*OX 
YL0SS=T*S+XPRIN2-X (21 
GO TO 62 

52 CONTINUE 

XL OSS =7 *R+XPRI Nl-X ( 11 
YLGSSM l.-FINALl*QY 
62 CONTINUE 

NOGO= 1 

X L 0 S ( J »K1=XL0SS 
Y LGS( J »K 1 = YLOSS 
RATL0S(J*K1 = Y LOS S -XLOS S 
FINTIMt J , K1=T 

IF( .NOT. (FATLOSU ,K 1 . LT.O -01 J GO TO 153 
SI GN=1.0 
GO TO 163 

153 CONTINUE 

SI GN--1.0 

163 CONTINUE 

R ATLOS ( J ,K1=ALAMDA*RATL0S ( J,K1 *( 1. 

€ - AL AM C A1 *S IGN*FIN TIM I J , K1 

1028 CONTINUE 

COUN T= CO UNT + 1 
I ED= 0 

GO TO 5000 
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6000 CONI INU E 

C 

C CALCULATE T1 

C 

IF ( • NOT .(TX.LT.TYiiGO TO 50 
T 1=T V-TX 
GO TC 60 

50 CONTINUE 

T 1=T X-TY 

60 CONTINUE 

C 

4 4 ***$ 4 X $ < $ 44c J$c 4>*£* $ ** $4 4<$* ❖*'*** * 4 4: 4c $ 4 ** 41 

C STAGE 2 

* $*4* *** * 4 * *•$*:££* 4>*******«*4t*****^^«4t^ 4* 4> $ 4 

c 

T ENC=0 . 0 
1 = 0. C 
CQUNT=1 
I ND= 1 

5 IFt.NOT.t (TEND.LE.T1J .AND. (NOGG. EQ.O )J ) 

£ GO TO 6 

T END =FLO AT (COUNT J*DEL 

IF( .NOT. (TX.LT.TYJ JGO TO 51 

CALL DVERK (N , FCN1 , T , X , T EN C ,TGL ,IND 
£ » CC tN W » W 1 1 ED J 

I F( .NOT. ( (X( 1 J.LE . (FINAL*CX)) .OR. 

€ ( X(2) .LE. (FINAL*QY) J )) GU TO 1C11 

xj n ^ = c x— x m ) 

YL0SS= (T*S J+XPR1N2-* (2 J*(TX<SJ 
RAT LOS ( J, KI=YLCSS-XLGSS 
XLOS( J ,KJ=XLOSS 
YL CS( J »Kj =YLOS S 
FI NTI M ( J , K J =TX +T 

IF ( .NOT •( RATLQS(J » K J .LT .0 .0 J ) GO TO 
& 15 A 



154 

164 

& 



1011 

51 

£ 

€ 



£ 



15 5 

16 5 

£ 

1013 



SI GN=1 . 0 
GO TO 164 
CONTINUE 
S I GN=- 1.0 
CONTINUE 

RATLO St J»K )=ALAMOA*RATLCS( J »K J+ 
(1 .-ALAMDA i*SI GN*F INTI M (J » K J 
NO G0= 1 
CONTINUE 
GO TO 61 
CONTINUE 

CALL DVERK (N » FCN2 » T , X , TEN D ,T0 L , I ND 
, CC tN W » W i I ED J 

IF( .NOT. ( (Xtli.LE . (FINAL* CXJ I .OR. 

( X( 2) .LE.(FINAL*QYI) )I GO TO 1013 
XLOSS=(T*R H-XPRiNl-XU J*(TY<RJ 
YLOSS=CY-X (2) 
RATLOS(J|KJ=YLOSS->LOSS 
XLQSt J»KJ = XLOSS 
YL QS( J »KJ = YLOS S 
FINTIMt J.K)=TY+T 
IF (.NOT. (RATLO St J»Ki .LT.0.0 JJGO 
TO 15 5 

SI GN=1 • 0 
GO TO 165 
CONTINUE 
SI GN=-1 .0 
CONTINUE 

RAT LOS ( J,K J = ALAMDA* RAT LOS (J,Ki + 
(1 .-ALAMDA J*SIGN*FINTIM(J»KJ 
NO G 0=1 
CONTI NUE 
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61 



CONT INUE 
COUNT-CUUNT+l 
I ED= 0 
GO TO 5 
6 CONTINUE 

C 

£44444 444444444 4 4 44c 4* 4 4444444444444444 4444444444444444444 

C STAGE 2 

£44444 444444444 4 4 4 4 44*4* 4444 4 4* 44444 4444 444 4** 4 44 4 44= #4= 4 44444 

c 



50 0 
1005 



£ 

£ 



1 = 0.0 
C QUNT= 1 
I ND= 1 
If (.NOT 
CONT 
T 
C 

t 

I 

( 



156 

166 

1014 



600 

C 

£44444 

c 

£*4c*44 

c 



I F( 
(X( 
GC 



CONTINUE 



. (NOGC.EG.O) J GO TO 600 
INUE 

END= FLOAT (COUNT i 
ALL CVERK ( N » FCN3 
CC »N W » W » I ED J 
F ( .NOT. ( ( X( 1 ) . L E 
X( 2) .LE. (FINAL*Q 
XLGSS = CX-X (1 ) 
YLCSS=CY-X(2) 

RA T LOS ( J » K ) =Y L 
XL 0S( J » Ki = XLOS 
YLOS( J,KJ= YLOS 
FI NTI M ( J . K ) =AM 
IF ( .NOT. ( RATLO 
TO 156 

SI GN=1 • 0 
GU TO 166 
CONTINUE 
S I GN=- 1 .0 
CONTINUE 
RA TLUS ( J » K ) =AL 
(1 .-AL AMUA )*SI 
NO G0= 1 
ONTI NUE 
OUNT = CGUNT+ 1 
ED=0 

NOT. ( (X( 1J .LE. ( 

) .LE • (FINAL*QYJ 
0 50 0 



^ CEL 

* 1 » X » TENC »TOL iIND 

. (F INAL* CX) J .CR. 
YJ ) J J GO TO 1014 



CSS-XLO SS 
<; 

S 

AXl ( TX » T Y I + T 

SI J, Ki.LT.0.0 JJGO 



F 

i 



AMDA4RAT L CS ( J ,K ) + 
GN*F 1NT Itf (J ,K j 



INAL*CX) ) -OR. 
JJ GO TO 1005 



♦4 44 4 44 44*444 444 44 44 444 44444444 4* 4*4=4: 4= 4- 4 44 4 44*4*$ 444444 

END OF eATTLE 

4444 444444444 44 4 44 44 4 4* 4 4 4*4 4444 444444444 44 44 444444 



XA ( 2CQUNTJ-XPR INI 
YA ( ZCGUNT) =XPR IN2 
NOGO s O 

ZA ( JC J -RAT LOS( J,KJ 
101 CONTINUE 

100 CONTINUE 

£44444 4444 4 4444 44444444444 4444 444444444444444444444444444 
C COMPUTE ROW MINIMUMAND SMALLEST kCWMIN 

£444444444 444 4 44444 44444 44444444444444444444444 4444444444 

C 



DC 301 KK=1,NS1 

RC WMI N ( KK J =R AT LOS ( KK# 1 ) 

DO 302 JJ=2, NS2 

If (.NOT . (RATLOS(KK, JJJ .LT.RClnMIN (KK i J) 
£ GO TO 1017 

RCWM IN (KKJ=RATLOS(KK , JJ J 
1017 CONTINUE 

302 CONTINUE 
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301 



CONT INUE 

DO 305 KM=2 ,NS1 

IF {.NOT # ( R GW MI N ( KK ) *LT #SMMIN)JGO TO 1019 
SMMIN=RCWM1 N (KK) 

C IN SER=K K 

1019 CONTINUE 

305 CONTINUE 

C 

Q*** ********** ************* ** ********* ******************* 

C SAMPLE 15*15 MATRICES FROM 60*60 MATRIX 

£***** **** **** ****** *************** ********************** 

c 

CG 30! JJ = i, 15 
J = 4* JJ 

DO 3C 8 KK= 1 » 15 
M4*K* 

RA T A ve ( JJ .KK l-RATLGS ( J , K i 
PAY( JJ | KK)-=RATAVE(JJ,KK) 

XLQSRR ( J J ,KK)=XLCS( J,K) 

UOSPKI JJ ,KK)=YLCS( J,K) 
f 7PR( JJ , KK) =F1 NT IM( J,K) 

308 CONTINUE 

307 CONTINUE 

C 

£****♦ *❖**>: *****^*** ******** *********** ************* ****** 

C OUTPUT MATRICES 

********* ************** * *** ************** ********** 

C 

WRITE! 1 ,580 ( (RATAVE( JJ,KK) ,KK»1,1S) , JJ»1, 15 J 
WRIT E ( I ,9814) 

WRITE ( 1 ,9811) ( (XLQSPR (JJ,KK) ,KK=1 ,15) , JJ=1 ,15) 
WRITE (1,9614) 

WRITE ( 1 ,9812) ( ( YLOSPR ( J J , KK ) , KK= 1 , 15) , JJ= 1 , 15 ) 
WRITE ( 1 ,9814) 

WRITE( 1,9813) ( (FT PR (JJ,KK), KK= 1,15), JJ=1, 15) 

WRITE ( 1 ,9814) 

980 FORMAT (' 1* ,15F7. 3//) 

9811 FORMAT ( ■ 1* ,15F7. 3//) 

9812 FORMAT (• 1* ,15F7. 3//) 

9813 FORMAT ( »1' ,15F7.3//J 

9814 FORMAT (•O* ,* • ) 

C 

Q********* **** ******** *********************************** 

C CALCULATE OP SOLUTION 

C ACD AES(SMMIN) TO ALL -PAY(J,K)TQ MAKE GAME VALUE 

C POSITIVE 

Q********* **************************** ******** * ********** 

c 

SMMI N = ABS ( SMM INJ+2. 

DO 403 J J= 1, 1 5 
RHS ( JJ)=1. 

08CCEF ( JJ) =1. 

00 4C4 KK* 1 , 15 

FAY( JJ, KK)= (PAY! JJ, KK1 + SMM1N ) 

404 CONTINUE 

40 3 CONTINUE 

C 

C CALL SUBROUTINE LP TO COMPUTE OPTIMIZED 

C STRATEGIES 

C 

CALL LP(PAY,RHS»QBCGEF,SMM1N) 

C 

c********* ********************* *********** ********* ****** 

C PLOTTING ROUTINE 

C***** ******* ********* **** *** ********************** ****** 

c 

CALL TEK61E 
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N FLQT= 1 

DO 102 L-l.NPLOT 

CALI FAGE( 8.50,11.0) 

CALL NCERDR 
CALL ELCWU P (1.0) 

CALL AREA2 D<5. 50,7.00) 

R EA C ( 5 , *) X VU 
WRITE (6,906)XVU 
R E A C ( 5 , *) Y VU 
WR ITE (6,907) YVU 
REAC (5,*)ZVU 
WR 1 1 E ( 6 ,90 8) ZVU 
C CALL FRAME 

CALL S CMP LX 
CALL X A XAN G ( 45 • 0 ) 

CALL YAXAN G( 45 .0) 

CALL ZA/AN G ( 45 • 0) 

CALL X3NAM E('X-CEPLOYMENT$' ,100) 

CALL Y3NAME( ' Y-CEPLOYM ENTS ' ,10C) 

CALL Z3NAMEI 'PAYOFF TO XS‘,100) 

C 

C INSERT LEADING IN ******* , REMCVE "C" 

C 

C CALL HEADI N ( ************** $ • , ICO, 1 • G , 4 ) 

C CALL HEADI N J •*************$ • , 100,1,0,4) 

C CALL HEADI N( '*************S' , 100,1.0,4) 

C CALL LEAD INC********* ****$' , 100, 1.0,4) 

C CALL MESSAGf*************!' ,100,1.2, 

C £ 7.0) 

CAi L VQLM3 C ( 1 . ,1. ,1.) 

CALL VUABS ( XVU ,YVL,ZVU ) 

CALL GRAF3D(2. 50, 'SCALE' , 7 .000 , 2.50 ,* SCALE* , 

C 5. 75, -5.0, 'SCALE* ,5.000) 

CALL BCX30 
CALL RASPL N ( 0. ) 

CALL SURMA T (KATLOS, 1.NS1, 1 ,NS2 ,0) 

CALL ENCPL(O) 

102 CONTINUE 

CALL OCNEPL 

906 FORMAT ( '0 * , 3X , 'XVU-' ) 

907 FORM AT (• O' ,3X , 'Y VU=' ) 

908 F0RMAT('C',3X,'ZVL=' ) 

C WRITE ( 1 ,988) I NSER , INSEC 

C988 FORMAT ( ' 0 ' , 3 X , ' I NSER= ' ,14,' INSEC=' ,14) 

WRITE ( 1 ,987) A ,B,C ,C 

987 FORMAT ('O' ,2X ,' A, E,C,C=' , 4F13 .3) 

WRITE ( 1 , 9 86) U , V , A L AMOA 

986 FORMAT ( • 0 • , 2X , • U , V - • ,2F 13.3, 2 X, ' L AMQA= 1 , F 1 3.2) 

C307 CONTINUE 
STGP 
END 
C 

C***** ******** * * ****************************** * ********** 

C SUER CL TINE FCN1 

C* ******** **** ******** **** ************ ********* ********** 

C 

SUEROUTINE FC N1 (N , T , X , XPR IME ) 

INTEGER N 

REAL X(N)»XPRIME( N),T 

COMMON /PAPAM/ A,C,U,V,S,R,B»D 

XPRIME(1)=-X( 1)*( U+A*X (2 ) )-B*X(2) 

X PRIME ( 2 A =-X (2)*( V*C*X(1) )+S-D*X( 1) 

RETURN 

END 

C 

c************* ******** **************** ******** *********** 
C SUEROUTINE FCN2 
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c# 444 4 4 444 4 44 4 4 4 4 44 44 4 4444 44^44444444444444444 $ $44 $ 444444 

c 

SUBROUTINE FC N2 (N ,T,X,XPRIMEJ 
INTEGEF N 

REAL X (N) ,XPRIME( N),T 

COMMON / P ARAM / A, C ,U , V ,S ,R ,B » 0 

X PRIME (1)=-X( 1I*( L*A*X (21 )+R-B*X(2J 

XFPIME (21 =-X(2)*( V*C*X 11 J )-D*X( II 

RETURN 

ENC 



C 

£4*44444444444 * 4 4 44 4:? 4 44 44 4 44 * 444 4 44444444 4 * 44 4 4444444 444 
C SUBROUTINE FCN3 

C*44 44 4 444 444444444 *4 44444444444444444444444444*4444444 44 

c 

SUBROUTINE FC N 3 (N , T, X , XPR IME J 
INTEGER N 

REAL >(NJ,XPRIME(NJ,T 

COMMON /PAPAM/ A. C ,U , V , S , R . B , D 

XPRIME(1)=-X( II *{ L+A *X(2) J-B4X(2) 

X PRIME (2J=-X( 2J*( \/+C*X(l) J-D*X(1J 
RETURN 
END 
C 

C**#*# 4444 4444 44444 44 444444 44444444444 44444444 4 4444 444444 

C SIBROLTINE FCN4 

£4 44 44 44444444 444*444444444444444444444444444444444 4 4 4 4 44 

c 

SUBROUTINE FC N4 (N , T, X , XPR IME J 
INTEGER N 

REAL X(N) ,XPR1ME( N),T 

COMMON /PAPAM/ A , C ,U * V , S » R , B , 0 

X PRIME ll)=-X( 1)*( U*A*X(2) I+R-B*X(2) 

X PRIME (2 J=-X(2)*( WC*X( 1J 1+S-C*X( II 

RETURN 

ENC 



£44444 4 444 444 4 4 4 44444444444444444444444444444444444444444 

C SUBROUTINE LP 

£4 44 44 44444444 4 4 4 4 4 4** 444* 4444444444444444444444444 444444 

c 

SUBROUTINE LP ( ALP , ELP ,CLP ,SMM INI 

INTEGEF IALP,NLP,M1,M2,IW (631 , I ELP 

REAL ALP (17,3 21 ,SLP,PSOL( 151 ,CSGL(15) ,RW(3E3J , 

£ SMMIN ,VAL ,BLP ( 15) ,CLP<15) 

N L F= 1 5 
M 1 = 1 5 
M 2=0 
I A LP= 1 7 

CALL ZXALPIALP, IALP , BLP,C LP, N LP ,M 1 ,M2 ,SLP, FSOL 
£ , CSOL ,FW , IN, I ELPI 

V A L= ( l./SLPI-SMMI N 
WRITE ( 1.725IVAL 

729 FORMAT (‘O' VALUE OF G AME= 1 , F 14 .3 I 
DO 405 J J = 1 , 1 5 

PSCL ( JJI=PSQL( JJI/SLP 
DSCL ( JJJ=D SOL( JJI / SLP 
WR 1 T E ( 1 ,7301 JJ ,FSCL(JJ I ,JJ ,CSOL( JJ I 

730 FCFMAT ( *0* , 2X, ' FRCB OF.STR. * 1 , 13 , 'OF V=' 

£ ,F 12.3, 'PROB.OF STR. #',13, 'OF X=',F12.3) 

405 CONTINUE 

RETURN 
END 
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APPENDIX I 



TRANSFORMING MIXED STRATEGY PROBLEM INTO LINEAR PROGRAMMING 



The payoff function has been defined to be 



A ( X , Y ) = L - L. 



Y selects his optimum mixed strategies which veild 



,, mm 
V = } max 



m m m 

T~fu q j , yi i 2j q i 23 



a mj q j 



J=1 



J = 1 



J = 1 



where 



a . . 
ij 



payoff to x when x adopts ith strategy and y adopts 
j th strategy 




probability that y selects j th strategy 



Subject to : 

j — 1 > Qj ^ d , j — 1 , 2 , . . . ,m 

j 

Let v Q ^ max S a 2 » • • • . 

Then the original problem becomes 
minimize Vq 
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Subject to : 



2Xi q j 

t'v'i 



< V 



< V 



0 



0 



E a . q . < v n 



= 1 



q j > 0 , j = 1,2, . . . ,m 



Dividing the constants by Vq (>0), we have 



1 m 
— > a, . q - < 1 

V 0 j7\ ^ 

i f^a,.q. < 1 
v n Z^ 2j { j 

M 



v o2 a mj q j 

U J* 1 m 
rn 

IBm ■ 



< l 



°» 3 



Let Q. = 



n a 

— and since 
v 0 



ram Vg = max — 



0 



= max [Q x + Q 2 + • • • + Q m J > 



the problem can be written as : 



max Q 0 = [Q 1 + Q 2 + • • • + Q m l 
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Subject to : 



E a i^j > 1 
j 

Z a 2jQ j > 1 

j 



Z a .Q. > 1 

j 



Qj > ^ » j 



= 1 , 2 , 



. , m 



1 4 i 

Since Q n = — and Q. = — 

0 v^ v 



0 



0 



= = = > qj - <y 0 - Q 



Q j 



0 



After solving the LP problem, the optimum strategies for y 

* rt 

is given by Vq Some constants, K = |min(a^j)| could 

have been added to a. . to ensure v n > 0. If this is done, 

K has to be subtracted from the optimum value obtained by the 
LP, that is 

v* = Q 0 - K 



where v = the value of the game. 
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APPENDIX J 



BACKGROUND DATA ON KOREAN WAR 

i: PERIOD CONSIDERED : 25 June 1950 to 7 July 1950 

(No American ground involvement 
as yet) 

2. TYPE OF ENGAGEMENT : Predominantly land combat 

3. GENERAL STATE OF READINESS : 

NORTH KOREA : Well prepared by 1950; arms 

build-up and training of 
troops since 1945; many 
military leaders and combat 
personnel were war veterans 
fighting in China 

REPUBLIC OF KOREA : By 1950; a small defense 

force began to take shape 
through American aid; 
training only started around 
1948. 

4. SOURCE OF DATA : a) Appleman, R.E., United States 

' Army in the Korean War , Department 
of the Army, 1961 . 

b) Montross, Lynn, U.S. Marine Opera - 
tions in Korea, U.S. Marine Corps., 
1954. 

5. RELATIVE STRENGTH : 

NORTH KOREA REPUBLIC OF KOREA 



a . 


Total strength = Q 

X 


135,000 men 


Qy = 95,000 men 


B. 


Tanks : 


150 


nil 




Artillery pieces : 


1,600 


700 


c . 


Aircraft 
(i) fighters 


- 40 


no combat aircraft 




(ii) attack bombers 


- 70 


(22 trainer, 4 




(iii) reconnaissance 


- 10 


auxiliary; no pilot) 
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6 . 



CORRESPONDING PARAMETERS USED IN MODEL 



a = 0.7 
b = 0.4 
u = 0.15 

1 unit = 13,500 men 



c = 1.0 
d = 0.6 
v = 0.2 
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